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BASIC  NOMENCLATURE 


Capital  letters  denote  matrices. 

Lower  case  letters  denote  scalers  and  constants. 

Underlining  denotes  vectors. 

Denotes  i**1  element. 

Denotes  a  nxm  matrix  of  all  zeros  with  a  1 
in  location  i,  j. 

Characteristic  equation  of  the  sensor. 

Low  frequency  2nd  order  term  in  A^. 

High  frequency  2nd  order  term  in  A^. 
Characteristic  equation  of  the  observer. 


ABSTRACT 


The  automatic  controls  division  of  the  Naval  Air 
Development  Center,  (NADC)  of  Warminster,  Pennsylvania,  is 
currently  investigating  the  use  of  the  Model  8160  rate  sensor 
offered  by  the  inertial  division  of  Systron  Donner,  Concord, 
California.  This  rate  sensor  appears  to  have  the  potential 
for  reduced  maintenance  cost  over  current  models,  however,  its 
output  response  has  zero  transmissibility  at  zero  frequency 
and  NADC  wanted  to  determine  if  this  could  be  altered  through 
use  of  a  filter. 

An  observer  was  designed  to  function  as  a  state  estimator 
for  the  rate  sensor.  Once  the  system  states  were  available, 
the  output  of  the  observer  was  modeled  after  a  desirable 
response.  The  resulting  observer  output  performance  exceeded 
expectations. 

This  report  is  a  sensitivity  study  of  the  effect  of  all 
the  rate  sensor  and  observer  parameters  on  the  generated 
observer  output. 
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1.0  INTRODUCTION 


1. 1  Background 

The  automatic  controls  division  of  the  Naval  Air  Develop¬ 
ment  Center,  (NADC)  of  Warminster,  Pennsylvania,  is  currently 
investigating  the  use  of  the  Model  8160  rate  sensor  offered 
by  the  inertial  division  of  Systron  Donner,  Concord,  California. 
Appendix  A  is  the  data  sheet  on  the  rate  sensor.  This  rate 
sensor  appears  to  have  the  potential  for  reduced  maintenance 
cost  over  current  models,  however,  its  output  response  does 
not  meet  the  requirements  of  NADC. 

In  a  previous  study,  [1]  an  observer  was  designed  to 
function  as  a  state  estimator  and  then  the  states  were  re¬ 
combined  so  as  to  achieve  a  desirable  output. 

This  report  is  a  continuation  of  the  previous  one.  The 
objective  is  to  investigate  the  sensitivity  of  the  generated 
output  as  a  function  of  the  rate  sensor  parameters  and  of  the 
observer  parameters. 
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2.0  SYSTEM  EQUATIONS 


Clearly  the  coefficients  b,  c,  d  and  e  are  functions  of 


^1*  wi»  ^2  an<^  (d2‘ 

The  transfer  function  of  equation  (2.1.4)  may  be  written 
in  phase  variable  canonical  form  [2]  as: 


0  10  0 
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q2 

q3 

Lq4 


or 


^  *  Ac^  +  B0 


eo  *  C3* 

The  q  states  are  phase  variable  states  and  are  defined 
as  *  <*2*  ^2  *  q3*  etc,»  t^ie  rate  sensor  output,  eQ  is  a 
scaler  signal. 


The  signal  6  is  not  available  as  an  input,  but  6  is. 


and  so  equation  2.1.4  may  be  written  in  terms  of  6  as 


>■ 


as 

~ Z - X - 2 - 

s^+es^+ds^+cs+b 


(2.1.6) 


This  transfer  function  has  a  set  of  states  different  from 
the  q  states.  Let  these  be  the  x  states  so  that  the  phase 
variable  canonical  state  equations  become: 
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2 


x 


3 
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4 


0  10  0 

0  0  10 

0  0  0  1 

-b  -c  -d  -e 

e  -  [0  0  0 
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(2.1.7) 


or 

*  »  Ax  +  B0 

e  *  C  .x 
o  4- 

Note  that  the  A  §  B  matrix  are  the  same  but  in  one  case  e 

o 

is  equal  to  the  qj  state  and  in  the  other  eQ  is  equal  to  the 
Xj  state.  This  is  a  consequence  of  the  phase  variable  definition 
of  states.  It  is  seen  that  the  relation  between  the  states 
is  i  ■  q. 
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2.2  Observer  Equations 

The  signal  b  is  not  available  and  so  the  observer  must 
use  8  as  its  input  reference,  hence  the  observer  is  defined 
[3,4,5]  as 


where  F  is  unknown, 
rate  sensor  output. 


2  -  F2  +  B6  +  Ke0  (2.2.1) 

B  is  the  same  as  in  2.1.7,  e„  is  the 

o 

and  K  is  the  vector  of  observer  gains 


If  the  rate  sensor  is  defined  by  equation  (2.1.5),  the 
error  between  the  sensor  states  and  observer  states  is  e  «* 
or 

e  ■  £  -  z  *  A<j  +  B0  -  Fz^  -  B0  -  Ke0  (2.2.2) 

A  problem  arises  in  that  0  and  §  can  not  be  compared;  hence 
it  is  necessary  to  use  equation  (2.1.17)  to  define  the  rate 
sensor  [1].  Using  (2.1.7),  e  •  x  -  z  and  so 

e  *  x  -  i  -  Ax  +  B0  -  Fz  -  B8  -  Ke0  (2.2.3) 

Using  eQ  from  (2*.  1.7)  gives 

4  *  k  -  4  *  Ax  -  KC4x  -  Fz  (2.2.4) 

If  F  *  A  -  &C4,  then 

e  -  (A  -  KC4)x  -  Fz 
or 

t  *  F£ 
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which  converges  if  the  eigenvalues  of  F  are  all  negative, 
and  the  K  vector  is  a  set  of  arbitrary  gains  used  to  ensure 
this. 

Hence,  the  observer  must  be  determined  in  terms  of  the 
x  states;  however,  once  the  observer  dynamics  are  defined 
it  makes  no  difference  which  states  are  used  from  thereon. 

The  auxiliary  output,  or  generated  output  [1]  is 


(2.2.5) 


-  [0  .87  4.04  13.33)z 

Therefore,  the  observer  matrix  equations  are: 


(2.2.6) 


e  *  [0  m2  mj  m^]z 


Figure  2.2.1  illustrates  the  above  relations.  As  shown, 
under  ideal  conditions,  the  signal  6  and  the  signal  6  are 
related  as  0  *  /0dt.  When  this  is  true,  then  the  analysis 
using  the  model  is  generally  simpler.  If  the  mathematical 
relation  between  6  and  0  fails,  then  the  relation  in  the 
top  of  Figure  2.2.1  must  be  used. 
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2.3  System  Equations 


If  the  dynamics  of  the  observer  are  augmented  to  those 
of  the  rate  sensor,  then  a  set  of  composite  states  may  be 
defined  as: 

v  -  Hv  +  L6  (2.3.1) 
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e  *  Dv 
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(2.3.2) 
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Equations  (2.3.2)  are  the  full  set  of  state  equations  that 
describe  the  dynamics  of  the  rate  sensor  and  observer. 


!  E 


The  solution  to  equation  (2.3.1)  is 
v  (s )  =  (sI-H) ' 1L6 (s) 


S(s)  =  D(sl-H)  ’  aL0  (s) 


(2.3.3) 


In  terms  of  equation  (2.3.2),  the  (sI-H)  matrix  is 


A  0 


kc4  f 


sI-A 


-KC4  sI-F 


(2.3.4) 


The  inverse  of  (2.3.4)  is 


(sI-H) 


(sI-A) 


l31'nj  "  |_(sI-F)  iKC4(sI-A)  x  (sI-F)'AJ 
The  A  matrix  was  defined  in  equation  (2.1.6)  as 


(2.3.5) 


0  10  0 

A  =  0  0  10 

0  0  0  1 

-b  -c  -d  -e 


and  so 


’32  2 

s  +es  +ds+c  s  +es+d  s+e 


(sI-A)'1  *  A^  'b 


3  2  2 

s+es  +ds  s  +es  s 

-cs-b  s^+es^  s' 

2  2 

-cs  -bs  -ds  -cs-b  s’ 


(2.3.6) 


where  AA  is  the  characteristic  equation  of  the  (sI-A) 


matrix  and  is 


A^  ■  s^+es^+ds^+cs+b 

2  2  2  2 
■  AA1AA2  ■  (s  +61u1s+ta)1)  (s  +262W2S+W2) 


# 


Next,  given  the  F  matrix  in  (2.2.3)  as 


0  -K. 


0  0  1  -K. 


0  1-K, 


(2.3.7) 


|_-b  -c  -d  -(e+K4)j 

then  equation  (2.3.8)  on  the  next  page  may  be  found. 

Knowing  the  relation  (2.3.5),  the  only  terms  missing 
in  (sI-H) * 1  are  the  (sI-F)'1XC4(sI-A)‘1  terms.  This  is  a 
4x4  matrix  and  in  general  the  terms  are  quite  complex. 
Because  it  is  encountered  quite  often,  it  is  defined  as 

Pll  C12  C13  C 1 4I 


(si  - F) "  aC^sI -  A)  * 


11 

C12 

C13 

C14 

21 

C22 

C23 

C24 

31 

C32 

C33 

C34 

41 

C42 

C43 

C44 

’  C34 

and 

C..  terms 
44 

;  are 

(2.3.9) 


often  and  so  are  written  out  for  future  reference: 


C24  =  K2s6  +  (K3+K2e)s5-*-(K4+K3e+K2d)s4-K1bs3 


C34  =  K3s6+(K3e+K4)s5+(-K1b-K2C)s4-K2bs3  (2.3.10) 

(^4  =  K4s6- (K1b+K2C-*-K3d)s5-(K2b+K3C)s4-K3bs3 
Additional  terms  will  be  found  as  needed. 
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+(e+K4)sJ+[d(l-K3)-cK2-bK1]st  +  [c(l-K3)-bK2]s  +  b(l- 


Referring  back  to  equations  (2.3.3)  and  using  (2.3.2) 


and  (2.3.5) 
written  as 

®(s)  = 


the  solution  of  the  composite  system  may  be 


(sI-A) 


[0  M] 


r 

Lsi- 


-i 


F)'1KC4(sI-A)'1  (sI-F)'1 


(2.3.11) 


£(s)  and  0(s)  are  scalers  and  so  the  remainder  of  the  terms 
are  the  transfer  function  of  the  observer  output  with  0  input. 
Expanding  (2.3.11)  gives 

g(s)  =  M{  (sI-F)  ~  1KC4  (sI-A)  "  1  +  (sI-F)  _1}B0  (s)  (2.3.12) 

The  expression  between  the  brackets  is  rather  messy, 
however : 

Theorem 

(sI-A) " 1  =  (sI-F) ”1KC4 (sI-A)  1  +  (s I -F) * 1  .  (2.3.13) 

Proof 

Post  multiply  both  sides  of  (2.3.13)  by  (sI-A)  to  get 
I  =  (sI-F)'1KC4  +  (si- F)‘1(sl-A) 
then  premultiply  both  sides  by  (sI-F)  to  get 
(si-  F)  =  KC4  +  (sI-A) 

But  by  definition,  equation  (2.2.4),  F  -  A-KC4  and  so 


sI-F  =  sI-A+KC4 
=  sI-(A-KC4) 
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Therefore,  equation  2.3.12  simplifies  to 


e  (s)  -  M(sI-A)'1B6(s)  (2.3.14) 

Of  course  this  is  not  valid  at  all  times  as  explained  below. 

2.4  Assumptions 

As  with  any  sensitivity  study,  the  object  is  to  determine 
the  effect  of  each  parameter  independently  of  all  others. 

Hence,  the  assumption  is  made  that  we  know  exactly  the  rate 
sensor  parameters  and  that  we  design  the  observer  with  exactly 
the  values  desired.  Then  the  effect  as  each  parameter  varies 
from  its  ideal  value  is  investigated. 

Along  with  this  is  the  assumption  that  8  =  /0dt,  the  signals 
0  and  0  are  assumed  to  have  the  proper  mathematical  relation. 

As  seen  with  the  Theorem  of  last  section,  if  the  0-0  relation 
holds  it  simplifies  the  analysis  at  times.  The  effect  of 
variations  between  0  and  0  is  also  investigated. 

In  some  cases  the  effect  on  the  output  changes  depending 
upon  where  the  parameter  changes.  That  is,  if  the  forward 
path  gain  of  the  rate  sensor  varies  10%  from  the  ideal,  the 
effect  is  different  than  if  the  forward  path  gain  of  the 
observer  is  off  10%. 

2. 5  Effect  of  the  Observer  Root  Locations 

The  characteristic  equation  of  the  observer  is  Ap  and 
is  defined  in  equation  (2.3.8),  where  it  is  seen  that  the 
roots  of  Ap  are  functions  of  the  observer  gains,  K.  The 
observer  is  appended  to  the  rate  sensor  and  the  total  system 


has  the  characteristic  equation  A^Ap,  where  A^  is  the 
characteristic  equation  of  the  sensor.  Therefore,  the 
observer  does  not  affect  the  roots  of  the  sensor,  but  adds 
its  own  roots  to  the  overall  system.  Since  the  roots  of 
Ap  depend  upon  the  observer  gains  K,  then  changing  the  gains 
changes  the  roots. 

In  this  case  the  observer  is  functioning  as  a  state 
estimator,  that  is,  it  is  generating  the  rate  sensor  states. 
The  observer  gains  are  arbitrary  to  the  extent  that  the  roots 
of  Ap  must  all  have  negative  real  parts.  After  this,  the 
magnitude  of  the  real  parts  of  the  roots  of  Ap  determine 
how  fast  the  observer  states  converge  and  track  the  sensor 
states.  As  found  in  [1],  placing  the  negative  real  part 
of  the  observer  roots  between  -2  and  -10  appears  to  give 
reasonable  results. 

From  the  standpoint  of  the  observer  tracking  the  sensor 
states,  the  observer  root  locations  do  not  matter  once  initial 
condition  transients  die  out.  That  is,  the  observer  root 
locations  determine  the  rate  of  convergence  to  the  plant 
states  initially,  but  they  all  track  the  same  from  then  on. 
Observer  roots  with  a  negative  real  part  of  -2  converge  to 
the  sensor  states  in  approximately  two  seconds,  and  track 
from  then  on.  Two  seconds  compared  to  the  start  up  time 
of  the  airplane  means  nothing. 


On  the  other  hand,  refer  to  Figure  2.2.1(a).  Note  that 
the  connection  between  the  system  and  the  observer  is  through 
the  gains  K.  If  the  gains  K  o,  then  the  output  8  is  totally 
dependent  on  the  observer  parameters,  which  reduces  to  the 
model  equations  of  the  sensor.  On  the  other  hand,  as  the 
gains  K  become  large,  dependence  upon  the  model  parameters 
decreases . 

In  general,  as  the  K  gains  increase,  the  real  part  of 
the  roots  of  the  observer  become  more  negative,  but  they  must 
do  this  according  to  Ap  of  equation  (2.3.8)  in  order  to  have 
a  stable  system.  To  insure  this,  the  roots  of  Ap  are  generally 
first  placed  and  the  K  gains  then  calculated. 

Throughout  the  study,  in  most  cases,  the  observer  gains 
had  effects  on  sensitivity.  In  all  cases,  it  was  assumed  that 
all  four  roots  of  Ap  were  at  -2,  -5,  or  -10  on  the  real  axis, 
and  thus  on  all  figures,  the  curves  are  identified  with  a 
2,  5  or  10.  The  corresponding  gains  are  [1]: 


Observer  Gains 


K1 

k2 

K3 

K4 

-2 

.  0253 

.0411 

.984 

-167.3 

Root  -  5 

-1.64 

|  2.34 

.39 

-155.3 

Locations 


-10 


-42.5 


41.3 


-8.77  -135.3 


3.0  SENSITIVITY  EQUATIONS 


3. 1  Background 

Sensitivity,  as  first  defined  by  Bode  and  as  defined 
in  most  basic  automatic  controls  texts  [6,7]  is  currently 
called  the  logrithmetic  sensitivity  ratio.  According  to 
this,  the  sensitivity  of  P(s)  with  respect  to  Q(s)  is 


„P  _  3P(s)/P(s)  31P(s) 

Q  '  3Q(s)/Q(s)  "  aTQCit 


or 


,P  -  Q(s)  .  3P(s) 
sq  Pliy  3QXiT 


Once  Sq  has  been  found,  in  terms  of  a  transfer  function,  then 
given  a  percentage  change  in  Q,  ^  :  ^  then  the  percent  change 
in  P  is 


3P 


~  AL 

p 


=  AS.s 
Q  bQ 


This  works  well  for  transfer  functions  of  relatively  low  order. 

P 

If  Sq  is  computed,  only,  as  a  time  solution,  note  that 

p 

the  solution,  P  is  in  the  denominator  of  Sq  .  In  addition, 
the  above  does  not  lend  itself  to  matrix  equations  in  the 
time  domain. 


3. 2  Logarithmic  Sensitivity  Trajectory  Function 

Therefore,  most  current  work  in  the  area  of  parameter 
sensitivity  studies  employ  the  sensitivity  trajectory  function 
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or  the  logarithmic  sensitivity  trajectory  function  [8-16] 


P  .  n3P  _  3P 


3Q7Q 


(3.2.1) 


P  • 


where  Q  is  a  parameter.  (Actually,  Wq  is  a  sensitivity 
coefficient  [8] . ) 

The  definition  of  sensitivity  as  given  by  equation 

(3.2.1)  is  the  one  used  in  this  report. 

In  equation  (3.2.1),  since 

9P  ~  AP 
3Q  '  AQ 

the  deviation  in  P,  AP  due  to  a  percent  change  in  parameter 
Q,  AQ/Q  is 

P  .  AP 
=  Q 


or 


iP  =  p 

Q  PQ 


(3.2.2) 


Thus,  given  a  percent  deviation  in  Q,  this  value  is  multiplied 

D 

times  the  value  of  to  find  the  deviation  in  P.  Note  that 
since  Q  is  a  parameter,  then  AQ/Q  is  in  general  just  a  number, 

p 

10%.  Therefore,  the  curve  Uq  is  actually  the  curve  of  deviation 
in  P;  only  the  magnitude  needs  to  be  altered.  Also,  note  that 
equation  (3.2.1)  is  easy  to  use  with  matrix  equations  and  time 
solutions. 


Equation  (3.2.1)  is  related  to  Sq  from 


and 


Sp  * 
bQ 


AP  =  P 


AP/P 

Tm 


% 
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so  that 


(3.2.3) 


3. 3  Sensitivity  Equations 

The  rate  sensor  dynamics  augmented  with  the  observer 
dynamics  were  written  as 


\ 


I 


i 


* 

a 


V  =  Hv  +  L0  (3.3.1) 

6  =  Dv 

The  solution  of  this  was 

Z(s)  =  (sI-H)'1L0(s)  (3.3.2) 

and  d(s)  =  D (si -H) ” 1L0  (s)  (3.3.3) 

Since  e  and  9  are  scalers,  then  (3.3.3)  may  be  written  as 

=  DCsI-H)'^  (3.3.4) 

and  the  right  side  of  the  equal  sign  is  the  transfer  function 
of  the  total  system  from  input  0  to  observer  output,  6.  Let 
the  elements  of  the  D  matrix  be  denoted  as  d^^ ,  and  h^  and 
1^  be  elements  of  the  H  and  L  matrices  respectively,  that  is 

11  -  thiji 


»  -  [diji 


Let  E1,J  be  a  matrix  of  order  nxm  where 

n,n  n,m 

every  entry  except  the  entry  at  row  i,  column 
unity.  Thus, 


has  zeros  at 
j ,  which  is 


2,1  = 

2, 


o 

l 


o  61 

o  o| 


In  the  case  that  E1  is  a  row  matrix  or  a  vector,  it 

n  ,m 

will  have  only  one  superscript. 

To  find  the  sensitivity  of  8  with  respect  to  an  element 
of  the  output  matrix,  D,  in  equation  (3.3.4),  it  is  necessary 


to  find 


S  3£  _  di  3e 

%  '  "*571  “53— 


(S.3.S) 


From  (3.3.4)  it  is  seen  that 


3  i 

~JcT 


-  =  E]  (sI-H)'1L0 
i  1,n 


and  so 


W*.  -  di  Ej>n  (sI-H)_1Le 


(3.3.6) 


(3.3.7) 


To  find  the  sensitivity  of  §  with  respect  to  an  element 
of  the  L  matrix,  first  from  3.3.4 


then 


inr  -  D(sl-»r14.ie 

3 


4.  •  *jiH*i-H)-1EJ>1e 


(3.3.8) 


(3.3.9) 


In  order  to  find  the  sensitivity  of  %  with  respect  to 
an  element  of  the  II  matrix,  first  from  3.3.4 


3e 

3K77 


D (s I-H)  ” 1  ^-(sI-H)‘1Le 


(3.3.10) 


then 


h. . D(sl -H) " 1  (sI-H)_1L0 

ij  anij 


(3.3.11) 
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Equations  (3.3.7),  (3.3.9)  and  (3.3.11)  are  each  expanded 
and  treated  independently  in  later  sections. 


4.0  TOTAL  SYSTEM  RESPONSE 


4.1  Frequency  and  Time  Responses 

The  combined  sensor  and  observer  state  equations  are 


v  *  Hv  +  L0 
fe  *  Dv 


with  the  s-domain  solution 

£(s)  =  DCsI-H)'1L0(s) 

Using  equations  (2.3.2)  and  (2.3.5)  gives 


(4.1.1) 


!{§}-[«  Ml  [(SI'A)'11  ,  0 

(sI-F)‘1KC4(sI-A)‘1  (sI-F)~‘ 


(4.1.2) 


*  M  {  (sI-F)‘1KC4(sI-A)'1  ♦  (sI-F)'1}  B  (4.1.3) 
Equation  (4.1.3)  may  be  reduced  by  use  of  the  theorem,  equation 


(2.3.13)  to: 


%  =>  M(sl  -A)  *  aB6 


(4.1.4) 


The  coefficients  of  M  were  deived  from  the  low  frequency 
second  order  terms  of  equation  (2.1.3)  [1],  and  so  equation 


(4.1.4)  reduces  to 


(4.1.5) 


s  +  26  2(*)2s+u>2 


This  was  the  response  designed  for  in  [1]  and  is  a  simple 

second  order  response.  Written  in  terms  of  6,  (4.1.5)  becomes 

2 


s  *26 2Ui2s*lli2i 


an 

e  (s) 


21 


(4.1.6) 


Figure  4.1.1  illustrates  the  frequency  domain  solution  of 
equation  (4.1.1)  with  a  sine  wave  input  at  6.  This  is  seen 
to  be  the  same  response  as  one  would  get  from  equation  (4.1.6). 

Figure  4.1.2  illustrates  the  time  response  of  the  output 
8  from  the  solution  of  the  state  equations  (2.3.1)  with  a 
unit  step  input  at  8.  Again,  this  is  the  response  one  would 
expect  from  equation  (4.1.6). 

If  the  initial  values  of  the  observer  and  the  sensor 
are  not  the  same,  then  there  is  a  period  of  initial  condition 
transients  that  must  decay.  Figure  4.1.3  illustrates  the 
response  of  £  due  to  different  initial  states  of  the  sensor 
and  observer.  The  input  6  was  a  square  wave  of  period  1  second 
and  unity  magnitude.  Once  the  initial  transients  decay,  the 
output  8  follows  the  input  S  according  to  equation  (4.1.6). 

The  rate  at  which  the  initial  transients  decay  depends  upon 
the  root  locations  of  the  observer. 

Figure  4.1.4  illustrates  the  same  solution  as  Figure 
4.1.3  except  that  the  initial  values  are  equal. 
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OF  OBSERVER  OUTPUT 


OUTPUT  RESPONSE  TO  UNIT  STEP  THETA  DOT 


FIGURE  4. 


TIHE  -  SEC 


5.0  OUTPUT  SENSITIVITY  WITH  RESPECT 
TO  THE  OUTPUT  MATRIX 


5.1  Sensitivity  Equations 

The  dynamics  of  the  rate  sensor  are  shown  in  Figure 

2.2.1  as 

eQ  =  C4(sI-A)_1Be  (5.1.1) 

and  the  observer  equations  are 


%  =  M(sl-F)*1(B0+Keo)  (5.1.2) 

where  M  is  the  output  matrix  and  is 

M  *  (0  m2  m3  m4]  (5.1.3) 

Physically,  M  is  the  summation  of  the  observer  states  and  was 
contrived  to  produce  the  desired  output,  thus  the  m^' s  are  the 
gains  of  a  summing  amplifier.  Since 


/\ 


from  equation  (5.1.2)  one  may  find  that 


(5.1.4) 


Be 

8m- 


I.  =  E11>4(sl-F)-1(B6+Ke0) 

Using  (5.1.1),  (5.1.5)  becomes 

||.-  Ej>4(sI-F)'1[I+KC4(sI-A)'1]B6 
which  with  the  use  of  the  theorem  in  (2.3.13)  becomes 

EU(SI-A)'lB6 


(5.1.5) 


(5.1.6) 


(5.1.7) 
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The  (sI-A)  and  B  matrix  are 

known  from  (2.3.6)  and 

(2.1.5) 

respectively, 

and  so 

equation 

(5.1.7)  becomes 

li  « 

a  E1 

’l  " 

3 

— 4 
A  1,4 

s 

s2 

e 

(5.1.8) 

S3 

Now  using  equation  (5.1. 

4)  gives 

8 

am_ 

% 

=  2 

6 

(5.1.9) 

£ 

sam_ 

y 

m3 

=  3 

aa 

0 

(5.1.10) 

€ 

2 

S 

=  s  am 
A 

is 

(5.1.11) 

aa 


where  is  the  characteristic  equation  of  the  rate  sensor. 

If  0  is  a  unit  step  input,  then  under  steady  state 
conditions 


€  _  am, 

m2 1  ss 


1 


and 


/V  A 


(5.1.12) 


Therefore,  deviations  in  the  summer  gain  m2  are  directly 
reflected  to  the  output  6.  If  e  is  a  unit  ramp,  implying 
acceleration,  then 
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ss 


1 

s 


(5.1.13) 


m. 


6 

|ss 


4.6 


0 


According  to  this,  at  steady  state,  with  6  a  ramp,  deviations 
in  8  are  4.6  times  the  percent  deviation  in  m^  and  the  deviation 
in  8  grows  with  respect  to  errors  in  m2.  It  would  seem  that 
this  condition  would  seldom  if  ever  exist  for  any  length  of  time. 

Figure  5.1.1  illustrates  the  frequency  domain  sensitivity 
of  8  to  the  m^  with  a  sine  wave  input  at  6.  The  peak  magnitudes 
are  essentially  the  same,  but  the  output  6  is  sensitivity  to 
each  of  the  itk  over  a  different  frequency  range.  As  expected, 
only  m2  has  an  effect  at  zero  frequency. 

The  curves  of  Figure  5.1.1  represent  deviations  in  §  for 

a  100%  deviation  in  the  parameters.  Thus,  if  m2  is  off  by 

+10%,  then  the  y^  curve  of  Figure  5.1.1  is  multiplied  by  0.1 
m2 

to  get  A8.  This  is  then  added  to  the  actual  6  response  of 
Figure  4.1.1.  Figures  5. 1.2,3,  and  4  illustrate  the  output 
response  e  to  +10%  deviations  in  m2,  m^  and  m^  respectively. 

The  way  the  y  curve  deviations  add  directly  to  the  8 
curve  is  seen  to  be,  ^actual  =  §ideal  +  Therefore, 

normally,  only  the  deviation  curves  will  be  presented 
throughout  the  report. 
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OUTPUT  SENSITIUITY  TO  THE  OUTPUT  MATRIX  ELEMENTS 


FREQ  -  RAD/SEC 


OUTPUT  RESPONSE  TO  +18*  ERROR  IN  HZ 


FIGURE 


OUTPUT  RESPONSE  WITH  +10*  ERROR  IN  M3 


FIGURE 


OUTPUT  RESPONSE 


Figure  5.1.5  illustrates  the  time  response  of  the 
sensitivity  equations  (5.1.9-11)  to  a  step  input  at  6.  The 
differences  in  the  frequency  ranges  are  also  apparent  in 
these  response  curves.  These  curves,  like  the  last  ones, 
illustrate  the  deviation  in  the  time  response  of  the  output, 

8,  to  1001  deviations  in  the  m^  parameters.  Thus,  if  m2 

a 

is  in  error  by  +10%,  then  the  deviation  in  §  is  A£  ■  O.ly 

in  2 

This  response  is  added  to  the  response  of  £  with  no  deviations. 

That  is,  ^actuai  s  ^ideal  +  ^*8ure  5.1.6  shows  the  output 
i  with  +10%  error  in  m2,  m3  and  m^,  one  at  a  time. 

Again,  if  one  takes  O.ly  from  the  sensitivity  curves  of 
Figure  5.1.5  and  adds  it  to  the  e  response  curve  of  Figure 
4.1.2,  the  result  is  Figure  5.1.6.  The  response  of  the  system 
to  the  step  input  at  ©  is  very  quick,  but  the  errors  in  the  mi 
coefficients  cause  long  term  transients. 


5.2  Results 


Figure  5.2.1  illustrates  a  solution  of  the  state  equations 
for  £.  Curve  1  is  a  response  with  +10%  error  in  m2,  which 
implies  that  e  should  go  to  a  steady  state  value  of  1.1,  since 

A 

^actual  =  (1+0,lwm  5®ideal.  Curve  2  rePresents  a  solution  of 

2 

8  under  ideal  conditions  except  that  0  is  a  unit  step  at  t  *  0, 
then  has  another  step  at  t  *  5  superimposed.  Curve  3  is  the 


same  as  curve  2  only  m2  is  in  error  by  +10%.  Note  that  errors 
in  m2  do  not  affect  the  tracking  performance  of  the  observer, 
once  transients  due  to  initial  start  up  have  decayed,  they 
do  alter  the  steady  state  value.  The  above  is  undet  the 


OUTPUT  SENSITIVITY  WITH  RESPECT  TO  THE  OUTPUT  MATRIX 


FIGURE 


OUTPUT  SENSITIVITY  WITH  RESPECT  TO  OUTPUT  MATRIX 


FIGURE 


! 


conditions  of  a  unit  step  at  6. 

If  0  is  a  ramp  then  the  situation  of  (5.1.13)  is  in 

existence.  Figure  5.2.2  illustrates  this  condition.  Curve 

* 

1  represents  the  ramp  input  6  -  0.2t  under  ideal  conditions. 
Curve  2  represents  the  same  input  with  m2  in  +10%  error. 
Clearly  the  error  grows.  Curve  3  represents  ideal  conditions 
with  0  =*  . 2t  until  t  ■  5  at  which  time  0  *  . 2t+l.  Curve  4 
represents  the  same  input  as  3  only  m2  is  in  error  by  +10%. 
Again,  the  error  grows. 


i 

i 
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OUTPUT  RESPONSE  TO  ERRORS  IN  M2,  THETA  DOT  »  STEP 


FIGURE 


TIHE  -  SEC 


6.0  OUTPUT  SENSITIVITY  WITH  RESPECT 
TO  FORWARD  PATH  GAIN 


6 . 1  Sensitivity  Equations 


Referring  to  Figure  2.2.1(b),  the  rate  sensor  equation 
is 

eQ  =  C4(sI-A)_1B0  (6.1.1) 

and  the  matrix  B  in  (6.1.1)  is 


where  a  is  the  forward  path  gain  of  the  rate  sensor.  Also 
from  Figure  2.2.1,  the  observer  equation  is 

%  -  M(sl-F)'1(B 8  +  Ke0)  (6.1.2) 

where  the  matrix  B  in  (6.1.2)  is  the  same  as  above,  but 
in  this  case  represents  the  forward  path  gain  of  the 
observer. 

Under  ideal  conditions,  these  would  be  the  same. 

However,  if  the  gain  of  the  rate  sensor  deviates  from  the 
ideal  value,  it  affects  8  in  a  different  manner  than  if 
the  gain  of  the  observer  deviates. 

The  logarithmic  sensitivity  trajectory  is 

4*aH  (6.1.3) 

in  either  case. 
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6.2  Sensitivity  to  Rate  Sensor  Gain 

To  find  /3a  of  the  rate  sensor  one  starts  with 
equation  (6.1.2)  and  finds 


A 


(6.2.1) 


and  then  from  (6.1.1) 


where 


3eo 
Ta~  " 


C4(sI-A)_1eJ  ^6 


Therefore  equation  (6.2.1)  becomes 


(6.2.2) 


(6.2.3) 


||  -  M(sI-F)'1KC4(sI-A)'1eJ>10  (6.2.4) 

If  this  is  multiplied  by  a  we  get  the  sensitivity  equation 

-  aM(sI-F)"1KC4(sI-A)'1Ej ^9  (6.2. S) 

Using  the  theorem  of  (2.3.13),  equation  (6.2.5)  may  be 
written  as 

uj  *  aM[  (si  -  A)  _1-  (sI-F)’1]  eJ^O  (6.2.6) 
The  (sI'A)’1  and  (sI-F)’1  matrices  are  given  by  (2.3.6) 

4 

and  (2.3.8)  respectively.  E4  ^  as  given  by  (6.2.3)  picks 
off  the  last  column  of  each  so  that  (6.2.6)  becomes 
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(6.2.7) 


y*  -  a(o;»2m^4] 


& 


r*  -1 

_1_ 

f14 

s 

2 

1 

f  24 

s 

3 

'  *F 

f34 

s 

1 

4k 

4k 

l - 

> 


where  the 


-1 


are  entries  in  the  (sI-F)  matrix. 


Equation  6.2.7  reduces  to 


AFAA2 


U4 


AF‘fl*4AA+AA2 


[  » 


2K2)  s-*m2lv3 


iTj« 


e  (6.2.8) 


Figure  6.2.1  illustrates  the  sensitivity  of  e  to  the 
rate  sensor  gain  in  the  frequency  domain  with  a  sine  wave 
input  at  0. 

The  output  is  quite  sensitive  to  the  sensor  gain, 
however,  the  sensitivity  decreases  as  the  observer  roots 
are  made  more  negative,  as  discussed  in  Section  2.5, 
implying  increased  K  gains. 

If  the  right  side  of  equation  (6.2.2)  is  expanded 
the  result  is 


s30 

*r 


(6.2.9) 


For  a  step  or  a  ramp  (implies  acceleration)  at  9,  the  steady 
state  value  of  3eQ/3a  *  o.  Hence  variations  in  the  rate 
sensor  gain  have  no  effect  on  the  steady  state  value  of  6. 
This  is  seen  to  be  the  case  in  Figure  6.2.1. 

Figure  6.2.2  illustrates  the  sensitivity  function 
in  the  time  domain  and  Figures  6. 2. 3-4  illustrate  the 
solution  of  %  with  ♦  101  deviation  in  the  rate  sensor  gain. 


OUTPUT  SENSITIVITY  TO  SENSOR  GAIN, 
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FIGURE  6.2.1  FREQ  -  RAD^'SEC 
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FIGURE 


OUTPUT  RESPONSE  TO  10V.  ERROR  IN  SENSOR  CAIN 


FIGURE 


OUTPUT  RESPONSE  TO  -10V.  ERROR  IN  SENSOR  GAIN 


FIGURE  6.2.4 


Additional  transients  are  introduced  due  to  the  mismatch 
of  gains. 

6.3  Sensitivity  to  Observer  Gain 

In  this  case,  the  gain  a  is  the  one  in  the  B  matrix 
of  (6.1.2);  hence  from  (6.1.2) 


and  so 


II  =  MCsI-D^eJ  j  6 


m|  -  a  MCsI-F)'1eJ_1S 


(6.3.1) 

(6.3.2) 


4  - 1 

Matrix  ^  picks  off  the  last  column  of  (sI-F)  from 

(2.3.8)  resulting  in 


or 


ya  * 


24 

E34 
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(6.3.3) 


=  ~  ^m3^3+zn2^2^s  *  ***2^3^  ^  (6.3.4) 


Figure  6.3.1  illustrates  the  frequency  domain  sensitivity 
of  the  output  due  to  deviations  in  the  observer  forward 
path  gain.  Note  that  there  is  a  steady  state  deviation 
in  S. 

In  equation  (6.3.3),  if  §  is  a  unit  step,  6  is  a 
unit  ramp  and  since 


*24 

-s2K2  -  s(K3-l) 

f34 

* 

-s2(K3-1) 

f44 

s3 

-  — 
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OUTPUT  SENSITIVITY  TO  OBSERVER  CAIN 


10“ 2  3  5  7j0-l  3  3  7  3  5  7 101  3  5  7 jg2 

FREQ  -  RAD/SEC 


F 

(■ 

t 


then  the  steady  state  value  of  is 


and  so 


Aft  -  — 
ss  a 


a1 


ss 


da  | 

a1 


ss 


(6.3.5) 


(6.3.6) 


Thus  if  the  observer  roots  are  at  -10,  and  if  the  observer 
gain  is  in  error  by  10),  the  deviation  is  0.1  of  the 
-10  curve  in  Figure  6.3.1.  That  is,  take  0.1  of  the  curve 
and  add  it  to  the  ideal  frequency  response  of  Figure  4.1.1. 

Figure  6.3.2  shows  the  time  solution  of  the  sensitivity 
equation  (6.3.4)  for  a  unit  step  input  at  §.  Figures 
6.3.3  and  6.3.4  illustrate  the  effect  on  the  response  of 
8  with  a  +  101  error  in  the  observer  gain  a. 

The  final  effect, or  the  steady  state  effect,  of  errors 
in  the  observer  gain  a  is  identical  to  the  effect  found  for 
output  matrix  parameter  m2  in  the  last  section.  This  is 
apparent  in  equation  (6.3.5). 

Figure  6.3.5  illustrates  that  the  error  in  observer 
gain  causes  some  additional  transients,  but  the  basic 
tracking  ability  still  exists. 
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FIGURE  6. 


TIME  -  SEC 


FIGURE  6.3.5  TIME  **  SEC 


7.0  OUTPUT  SENSITIVITY  WITH  RESPECT 
TO  THE  OBSERVER  GAINS 


7.1  Sensitivity  Equation 


The  rate  sensor  dynamics  are  defined  by 
e0  -  C4(sI-A)_1B0 

and  the  observer  dynamics  are  defined  by 


6  -  M(sl-F)‘A(B0+Keo) 


and  so  the  logarithmic  sensitivity  function  desired  is 
e  _  „  de 


(7.1.1) 


(7.1.2) 


„J.  -  (7.1.3) 

From  (7.1.2)  one  may  write 

If.  -  M(SI-F)‘1Bj1e0  (7.1.4) 

Combining  (7.1.1)  and  (7.1.4)  gives 

||  -  M(sI-F)'1Ej>1  C4(sl-A)"1B0  (7.1.S) 

Referring  to  equations  (2.1.7)  and  (2.3.6),  by  inspection 

- 1  3 

it  is  seen  that  the  term  C4(sl-A)  B  reduces  to  as  /AA,  so 


If.  -  M(SI-F)-1Bj1e0 


ff.  -  M(sI-F)'1Ej1  C4(sl-A)*1B0 


that 


38  as* 


M(sI-F)-1Ej>1S 


Equation  2.3.8  is  the  expression  for  (sI-F) 
this  as 


(sI-F) 


(7.1.6) 

Represent 


11 

f12 

f13 

21 

f22 

f23 

31 

f32 

f33 

41 

f42 

f43 

(7.1.7) 


and  since  a  *  [0  a  2 


as 


dd 

TK, 


as 


and  then  of  course 


M* 

*i 


aK.  s 


m3m4],  equation  (7.1.6)  may  be  written 

^jfji)S  (7.1.8) 

j-2 

i-1 

2  4 

_  (7.1.9) 

j-2 

i-1 


In  all  cases,  even  when  5  is  a  ramp,  the  steady  state 
sensitivity  is  zero. 

Figures  7. 1.1-4  illustrate  the  sensitivity  of  e  to 
the  various  observer  gains,  K.  In  general  the  system  is 
more  sensitive  to  gain  Kj  because  it  is  most  closely  tied 
to  the  output  of  the  sensor.  The  deviation  curves  of 
Figures  7. 1.1-4  may  be  added  directly  to  the  ideal  response 
curve  of  Figure  4.1.1.  Figures  7.1.5  and  7.1.6  illustrate 
this  where  the  observer  roots  were  all  placed  at  -10  and 
then  the  observer  gains,  K  were  deviated  by  ♦  10%,  one  at  a 
time. 

Figures  7.1.7-10  illustrate  the  time  solution  of  the 

i»v  equations  in  (7.1.9)  to  a  unit  step  at  §.  Figures 
Ki 

7.1.11-14  illustrate  the  output  response,  8  with  a  10% 
error  in  the  observer  K  gains.  Figure  4.1.2  shows  the 
response  with  no  deviations. 

Errors  in  the  observer  gains  alter  the  observer  output 
in  the  frequency  band  of  approximately  0.1-100  rad/sec, 
depending  upon  which  gain  is  considered.  The  errors  cause 
no  steady  state  deviation  of  8. 
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8.0  OUTPUT  SENSITIVITY  WITH  RESPECT 
TO  THE  INPUTS 


8.1  Relations  of  the  Inputs 

Referring  to  Figure  2.2.1(a)  the  observer  has  as  an 
input  the  signal  6  because  8  is  unmeasurable;  in  fact. 


the  function  of  the  sensor  is  to  generate  a  signal  equal 
to  8. 


In  order  to  study  the  effects  of  the  other  parameters 
on  the  output  response,  it  has  been  assumed  that  9  ■  /§ dt. 
In  actual  fact,  the  signal  8  is  instrument  generated  and 
so  may  deviate  from  the  ideal  mathematical  relation. 


8.2  Sensitivity  with  Respect  to  6 


From  Figure  2.2.1(a)  it  is  seen  that  8  only  enters 


into  the  observer  equation  as 


and  so 


8  -  M(sI-F)’a(B8  ♦  Ke  ) 


||  -  M(sl-F) -1B 


(8.2.1) 


from  which 

y|  -  M(sl-F) _1B8 

M  and  B  are  known  from  equations  (2.1.7)  and  (2.7.5)  and 
hence  rather  simply  one  finds 

v6  "  G»2f24  +  ®3f34  +  »4f44)0 


(8.2.2) 


(8.2.3) 
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which  reduces  to 

uf  -  ^  [tt4AAl  *  fr>3K3*m2K2)s  *  ®2K3^  (8.2.4) 

F 

Comparing  these  equations  to  (6.3.3)  and  6.3.4),  it  is  seen 
that  the  sensitivity  due  to  variations  in  6  are  identical 
to  variations  in  the  observer  forward  path  gain.  If  one 
thinks  in  terns  of  a  transfer  function  from  8  to  the 
observer  output  and  that  a  is  the  transfer  function  gain, 
then  this  result  is  obvious. 

8. 3  Results 

The  results  of  section  6.3  apply  directly  to  this 
case.  In  particular  note  Figures  6.3.1  and  6.3.4.  These 
represent  the  response  of  6  with  a  101  error  in  either 
observer  gain  a  or  6.  This  is  a  constant,  fixed  gain  or 
offset.  The  observer  output,  0  is  off,  and  there  are  some 
additional  transients,  but  the  basic  tracking  ability  still 
exists. 

A  problem  that  could  arise,  and  must  be  considered 
is  the  effect  of  noise  on  8  .  Figure  6.3.1  indicates  that 
the  output  is  very  sensitive  to  this  noise  in  the  frequency 
ranges  of  0.1  to  100  rad/sec. 

Figure  8.3.1  is  a  solution  of  the  state  equations  for 
0.  The  observer  roots  are  all  at  -10  and  the  input  0  is 
taken  to  be 

8  ■  /  §dt  ♦  0.01sin(10t)  (8.4.1) 


This  means  that  8  is  the  true  integral  of  §  but  has  some 
10  rad/sec  noise  superimposed.  The  effect  of  this  small 
magnitude  noise  is  apparent. 

Figure  8.3.2  illustrates  8  under  the  same  conditions 
except  that  6  has  a  random  signal  of  magnitude  .01 
superimposed. 
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9.0  OUTPUT  SENSITIVITY  WITH  RESPECT 
TO  PARAMETERS  OF  THE  SENSOR 


9. 1  Sensitivity  Equations 

From  Figure  2.2.1(b)  the  rate  sensor  equation  is 

eQ  -  C4(sI-A)'1B0  (9.1.1) 

and  the  observer  output  is  given  as 

0  -  M(sl-F) ’1(B0  +  KeQ)  (9.1.2) 

The  parameters  of  concern  are  those  in  the  A  matrix  of 
(9.1.1).  There  are  four  parameters,  6^,  62  and 

Denote  these  as  the  general  parameter  p,  then  by  the 
chain  rule 

38  38  3eo  3A 

TJp  7fe~  3  A  "Sp  (9.1.3) 


Application  of  (9.1.3)  gives 

38  .  ufc t  ft'1  k  3e° 

Jp  *  M(sl-F)  K  -gp- 

|§  -  M(sl-F) -1  K  C4  (sI-A)’1  ||  (sI-A) "*B0 


(9.1.4) 


The  A  matrix  is 


0 

0 

0 

-b 


1 

0 

0 

-c 


0 

1 

0 

-d 


0 

0 

1 

-e 


Where  b,  c,  d,  and  e  are  functions  of  the  parameters. 
Therefore, 


3A 

3P 


0 

0 

0 


0 

0 

0 


0 

0 

0 


0 

0 

0 


-3b  -3c  -3d  -3e 

Up  3p  ‘Sp  'SpJ 
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The  B  matrix  of  (9.1.4)  is 


and  the  (sI-A)~x  matrix  is  given  in  equation  (2.3.6). 

The  product  is  pL  ~ 

(sI-A) ”1B  -  %  a 

S. 


and  hence  the  product  of 


%  w-A)'1*  -  [oj^ 


,  r3b  .  3c  .  2  3d  .  3  3e> 

r  (-5p  +  s-5p  +  8  *p  +  s  *p> 


Since  the  may  be  removed,  (9.1.8)  becomes 


r  -  -  (b  ♦  sc  ♦  s2d  ♦  s3e) 


The  characteristic  equation  of  the  sensor  is 

■  s*  ♦  s3e  ♦  s2d  +  sc  ♦  b 
and  so  (9.1.9)  may  be  written  as 

r  "  '  -sf  AA 

and  (9.1.7)  will  then  become 

(si-A)  _ib  -  r  s  ~i  a 


3,  ^ 

Ip  AA 


(9.1.6) 


(9.1.7) 


(9.1.8) 


(9.1.9) 


(9.1.10) 


Next, 


so  that 


C4(sI-A)  ■  2“  I*4x  *42  *43 

V*!-*)'1  j£  OI-A)"1®  •  *4  ^  C-4a) 


The  (sI-F)’*  term  has  the  form 
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f12 

f13 

f14 
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f22 

f23 

f24 

31 

f32 

f33 

f34 

41 

f42 

f«3 

f44 

(9.1.11) 


and  since  «i-  l°m  2  ^  m43  then 

rv 

m  (sI-F)_1K  -  ^  D^f -2  E&^fjj  E^f^j  *2  (9.1.12) 


-1  1  4  4 

(sI-F)  AK  -  -4-  Z  (K.  E  m  .f.  .) 

*F  2-j  J  i-2  1 


lk4J  i-2-4 
(9.1.13) 


Conbining  equations  (9.1.11),  (9.1.13)  and  the  fact  that 

6 

Up  *  p36/3p  produces  the  desired  result 

4 -  (p!pV(#>  i }i«j  j2 (9.1.14) 


3 

The  tens  P  -jp  AA  depends  upon  which  of  the  paraneters  is 
considered. 
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The  characteristic  equation  of  the  rate  sensor  was 

given  in  equation  (2.1.13)  as 

2  2  2  2 
Aa  -  AA]AA2  "  (s  +262w2s+w2^ 

Therefore,  for  6^; 

P  Ip  4A  ■  4lt2“lsI  4A2  (9.1. IS) 

for 

P-jp  ^  ^a2  (9.1.16) 

for  . 

P  Va  "  62t2w2sl  AA1  (9.1.17) 

and  for 

p  Ip  AA  “  w2*262s  +  2w2^  AA1  (9.1.18) 

Figures  9. 1.1-4  illustrate  the  sensitivity  of  6  with 
respect  to  6^  u>, ,  6.,  and  u2  respectively.  As  expected, 
parameter  sensitivity  decreases  as  observer  gains  K  increase. 
The  output  6  is  sensitivity  to  c^,  where  u>2  is  the  high 
corner  frequency  of  the  rate  sensor  itself.  As  seen  from 
equation  (9.1.4)  none  of  the  parameters  have  any  steady 
state  effect  even  when  i  is  a  ramp. 

Figures  9. 1.5-8  illustrate  the  time  solution  of  equation 
(9.1.14)  for  6j,  Mj,  $2  and  resPectively*  As  is 
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OUTPUT  SENSITIVITY  TO  SENSOR  PARAMETER  OMEGA  2 


FIGURE  9.1 


generally  always  the  case,  the  transient  sensitivity 


increases  as  the  observer  gains  K  increase.  However, 

A 

there  is  no  steady  state  deviation  in  e  due  to  parameter 
error. 

A 

Figure  9.1.9,  10  and  11  illustrate  the  output  of  e, 
(observer  roots  at  -10)  for  6  equal  to  a  step  of  magnitude  two, 
which  drops  to  a  magnitude  of  one  at  t  *  5  sec.  Transients  at 
the  steps  increase,  but  again,  the  basic  tracking  ability 
remains . 
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10.0  OUTPUT  SENSITIVITY  WITH  RESPECT 
TO  PARAMETERS  OF  THE  OBSERVER 


10.1  Sensitivity  Equations 

From  Figure  2.2.1(b)  the  rate  sensor  equation  is 

eQ  -  C4(sI-A)'1B6 
and  the  observer  output  is 

%  -  M(sl-F)_1(B0  ♦  Ke0) 

The  parameters  of  concern  are  those  of  the  F  matrix  in 
(10.1.2).  These  are  the  four  observer  parameters  6^, 
Wj,  $2  and  u>2  •  As  in  the  last  section,  let  p  represent 
a  general  parameter.  Then  from  (10.1.2) 

ff  -  M(sl-F)"1  ||  (sI-F) "*(B0  ♦  Ke0) 

Using  equation  (10.1.1)  the  terms 

(sI-F)’1(B6  ♦  KeQ)  -  [(sI-F)_1+  (sI-F) "1KC4(sI-A) _1]B0 

-  (si -A) " 1B0 

by  use  of  the  theorem  (2.3.13).  Hence  (10.1.3)  becomes 
||  -  M(sl-F)"1  |E 
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(10.1.1) 

(10.1.2) 

(10.1.3) 

(10.1.4) 

(10.1.5) 


4 


t 


The  F  matrix  is 


0 

0 

0 

-b 


1 

0 

0 

•c 


0 

1 

0 

■d 


*K1 

-*2 

i-k3 

-(e*K4) 


and  since  b,  c,  d,  and  e  are  functions  of  the  parameters. 


3F 


0  0  0  0 

0  0  0  0 

0  0  0  0 

3b  3c  _  3d  _  3e 

T5j>  '  3p  ~3p  '  Tp 


(10.1.6) 


Comparing  (10.15)  to  (9.1.5)  it  is  seen  that  3F/3p  -  3A/dp, 
consequently,  from  equation  (9.1.10),  one  may  write 


||CsI-AJ'1B  - 


0 
0 
0 
-  3 


a 


(10.1.7) 


L  *p  A^| 

In  a  like  manner,  from  equation  (9.1.12)  one  finds  that 


M(sl-F)'1  -  jp  [EJhjfjj  E»i^i2  zmifi3  rBifi4l  (10.1.8) 

F  i  -  2-4 

Combining  equations  (10.1.7)  and  (10.1.8),  (10.1.5)  becomes 


If  ■  zp^  (  •  fp  V  £2  mifi«)e 


(10.1.9) 


A 

and  since  Up  »  p38/3p,  we  get 

(10.1.10) 

Note  that  (10.1.10)  has  8  as  its  input.  The  term 
p3A^/3p  is  identically  the  same  as  found  in  the  last 
section.  The  results  are  repeated  below. 

For  6^; 

plp  AA  *  261  “  1  sAA2 

for  (Djj 

AA  *  C2fii  “i  s  +  2wi  )  aa2 

for  62I 

*  2^2(02SAA1  (10.1.13) 

for  o>2; 

p^t  AA  “  (262“2  s  +  2w2)  AA1  (10.1.14) 

Figures  10.1.1-4  illustrate  the  sensitivity  of  6  to 
the  observer  parameters  «1»  «jt  *2  and  u>2  respectively. 

The  summation  term  on  the  right  side  of  (10.1.10) 
was  expanded  in  equation  (6.3.4).  The*  term  was  shown 
to  reduce  to  unity  under  steady  state  conditions  in 
equation  (6.3.5),  assuming  a  unit  step  input  at  §. 


(10.1.11) 

(10.1.12) 
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Therefore,  assuming  a  unit  step  input  at  0,  under  steady 
state  conditions,  equation  (10.1.15)  becomes 


ss 


lim 

S-+-0 


When  considering  either  of  the  6  terms,  equation  (10.1.11) 
or  (10.1.13)  applies,  and  Ug|ss  a  0*  This  is  shown  in 
Figures  10.1.1  and  10.1.3. 

When  considering  either  of  the  to  terms,  equation 
(10.1.12)  or  (10.1.14)  applies.  Using  (10.1.12)  gives 


j 

0). 


ss 


,  2  2 
2<D^a>2 

'"2"”  2 
“l  w2 


=  2 


which  also  holds  for  u>2.  This  is  seen  in  Figure  10.1.2 
and  10.1.3. 

Figures  10.1.5-8  illustrate  the  time  response  of 
the  parameter  sensitivity  equation  (10.1.10)  for  6^ 

62  and  w2  respectively. 

Figure  10.1.9  and  10  show  the  output  response  for  +10% 
errors  in  62  and  respectively,  (observer  roots  at  *10).  For 
comparison,  the  response  of  the  ideal  system  is  also  plotted. 
The  input  ,  0,  was  aunit  step  until  t  »  3,  at  which  time  it 
went  to  a  magnitude  of  2,  and  then  returned  to  unity  at 
t  ■  5.  In  both  cases  the  basic  tracking  ability  remains, 
however,  the  parameter  errors  cause  some  transients. 
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FIGURE  10. 


Figure  10.1.11  is  a  run  of  $  with  +S%  error  in 
and  is  a  very  poor  response.  Additional  runs  for  <5^  were 
unstable  for  even  +51  error  in  6^.  On  the  other  hand. 
Figures  10.1J2and  13  were  runs  made  with  -101  error  in 
oij  and  respectively.  Since  these  are  low  frequency 
terms,  the  transients  don’t  decay  rapidly,  but  the  response 
is  well  behaved. 

If  the  input  9  is  a  unit  ramp,  then 


lim  2 
s-+o  s 


oo 


This  means  that  the  error  grows  under  conditions  of  long 
term  constant  acceleration.  This  is  the  same  result  as 
discussed  in  Section  5. 

In  a  similar  manner,  if  0  is  a  unit  ramp,  then 


and 


which  means  a  small,  constant  error  under  constant 


acceleration. 
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11.0  SUMMARY 


In  general  a  sensitivity  study  involves  a  large 
amount  of  data,  and  this  one  is  no  different.  A  conden¬ 
sation  based  on  effects  is  attempted  below. 

When  the  velocity,  9  suffers  a  step  change,  errors 
in  the  observer  parameters  and  gain  a;  the  output 

matrix  element  »2*>  an<*  steady,  offset  errors  in  6  the 
generated  signal,  affect  the  steady  state  output  of  the 
observer, 

When  the  velocity  0  is  a  ramp,  constant  errors  in 
the  above  mentioned  parameters  tends  to  cause  ever  in¬ 
creasing  errors  in  the  output  e.  In  addition,  the 
output  matrix  element  and  observer  parameters  6^  and 
62  affect  the  steady  state  output  of  the  observer.  These 
would  appear  to  be  "second  order"  effects  and  will  most 
likely  not  appear  in  the  actual  aircraft. 

After  initial  condition  transients  have  died  and 
the  observer  has  reached  steady  state,  the  following 
parameters  affect,  to  some  extent,  the  performance  of 
the  observer  while  tracking  step  changes  in  0. 

Observer  gain  a 
Sensor  parameters  62  and  o>2 
Observer  parameters  6^,  u>j,  62  and  uj2 
Noise  on  the  generated  signal  0 


The  major  problem  appears  to  be  with  the  sensor 
parameters  Uj  and  6^.  If  these  are  in  error  on  the  low 
side  an  acceptable  response  results.  If  they  are  in 
error  on  the  high  side,  the  solution  is  unstable. 


12.0  REFERENCES 


1.  Mitchell,  E.  E.,  "Observer  Design  for  Improved  Rate 
Sensor  Performance",  ADB026070L. 

2.  Kuo,  B.  C. ,  Automatic  Control  Systems,  3rd  ed. 
Prentice-Hall,  1975. 

3.  Luenberger,  D.  G.,  "An  Introduction  to  Observers", 

IEEE  Trans  on  Auto  Control,  Vol  AC-16,  No.  6,  Dec.  1971. 

4.  Meditch,  J.  S.  and  G.  H  Hos tetter,  "Observers  For 
Systems  With  Unknown  and  Inaccessible  Inputs". 

Int.  J.  Control,  1974,  Vol.  19,  No.  3,  pg.  423. 

5.  Hostetter,  G.  H.  and  J.  S.  Meditch,  "On  the 
Generalization  of  Observers  to  Systems  with  Unmeasurable, 
Unknown  Inputs",  Automatica.  Vol.  9,  pg.  721,  1973. 

6.  Kokotovic,  P.  V.  and  Rutman,  R.  S. ,  "Sensitivity 

of  Automatic  Control  Systems  (Survey)".  Automation 
and  Remote  Control,  Vol.  26,  No.  4,  1965,  pg.  TTT7~ 

7.  Dorf,  R.  C. ,  Modern  Control  Systems,  Addison-Wesley , 

1967. 

8.  Brewer,  J.  W.,  Control  Systems  Analysis,  Design 
and  Simulation,  Prentice  Hall,  1974. 

9.  Cruz,  J.  B.,  Feedback  Systems,  McGraw-Hill,  1972. 

10.  Frank,  P.  M. ,  Introduction  to  System  Sensitivity, 

Academic  Press,  197fi. 

11.  Proceedings  of  the  3rd  IFAC  Symposium,  June  18-23, 

1973,  Ischia,  Italy. 

12.  Cruz,  J.  B. ,  Ed.,  System  Sensitivity  Analysis,  London: 
Hutchinson,  1973. 

13.  Ngo,  N.  T.,  "Sensitivity  of  Automatic  Control  Systems", 
Automatic  and  Remote  Control,  Vol.  32,  pg.  735,  1971. 

14.  Porter,  W.  A.,  "The  Interrelationship  Between  Observers 
and  System  Sensitivity",  IEEE  Trans,  on  Auto.  Control, 
Feb.  1977. 


L 

-* 


in 


15.  Bongiorno,  J.  J.,  "On  the  Design  of  Observers  for 
Insensitivity  to  Plant  Parameter  Variations",  Intr. 

J.  Control,  Vol.  18,  No.  3,  1973,  pg.  597. 

16.  Cooper,  D.  J.,  and  A.  Graham,  "Errors  in  State  Variable 
Reconstruction  When  an  Observer  is  Subject  to  Constant 
Input  Disturbances  and  Measurement  Noise",  IEEE  Trans. 
on  Auto.  Control,  Feb.  1977. 

17.  Mitchell,  E.E.  and  C.  F.  Olsen,  "Digital  Computation 
of  Frequency  Response",  Simulation,  Nov.  1975,  pg.  137. 


112 


INPUT 


APPENDIX  A 


SYSTRON-  CONNER 
MODEL  8 160  RATE  SENSOR 
TRANSFER  FUNCTION 


P/N  8160-120- PI5 


0.075 


6 


S4 


S‘+0.303S*I- 0.0656 


15600 _ 

1755  +• 15600/ 


CORRESPONDING  AMPLITUDE  RESPONSE: 


g-2 -T6 


INtRTIAL  OIVtaiON 

SVSTIIOM  OOMMBI 


APPENDIX  B.  COMPUTER  PROGRAMS  FOR  FREQUENCY 

RESPONSE  OF  SENSITIVITY  EQUATIONS 

A  sequence  of  main  programs,  all  essentially  the  same, 
were  written  to  generate  the  polynomials.  The  main  programs 
then  called  the  subroutine  SENSUB  to  generate  the  frequency 
response  curves.  The  algorithm  used  in  SENSUB  is  explained 
in  reference  [17] . 

The  programs  are  as  follows: 

A 

SENSE  -  computes  the  frequency  response  of  e. 

A 

Compute  the  frequency  response  of  e  with  respect  to: 
SENSEM  -  the  nu  of  the  output  matrix 
SENSEAO  -  the  observer  gain 
SENSEAS  -  the  rate  sensor  gain 
SENSEZ1  -  Sj  of  rate  sensor 
SENSE"^  -  W  of  rate  sensor 
SENSEZ2  -  S2  rate  sensor 
SENSEK  -  observer  gains 
SENSEW10  -  W  of  observer 
SENSEZ10  -  ^  of  observer 
SENSEW20  -  2  of  observer 

SENSEZ20  -  62  observer 
SENSUB  -  frequency  response  subroutine 
POLYPAK  -  polynomial  subroutines 

The  programs  SENSEZ1,  SENSEwl,  SENSEZ2  and  SENSEu2  are 
the  same  except  for  the  polynomial  forms.  Hence,  SENSEZ1 
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SENSE 


100*** 

120*** 

140*** 

160*** 

180  COtMON  GN,GD,IGN,IGD,WN,WI,IGNP,IGDP,CV 

200  IMPLICIT  DOUBLE  PRECISION ( A-H ,0*Z) 

220  DOUBLE  PRECISION  K1 ,K2 ,K3»K4,M 

240  INTEGER  IGN(20) ,IGD<20) ,WN 

260  DIMENSION  GN(2O,1O),GD(2Q,10) 

265  DIMENSION  F(3,4),H<3,7> 

266  DIMENSION  M(3),DELA1(3),DELA2(3),DELF(7) 

280  CHARACTER  CV(2,40) 

281  CHARACTER* 1  P1(40) ,P2(40) 

285  DATA  P2/3*1HA,1H*,1H*,1H8/ 

287  DATA  P1/3*"A",  , 2*"A" , 2*"*" , "S" , 3*"A" , 2**+* , 

288  &  "A",  "*«,"R","-t.","A", 

300  LIBRARY  "SENSUB" 

320  B2(4)*"*" 

340* 


MAIN  PROGRAM  TO  COMPUTE  THE  FREQUENCY 
RESPONSE  OF  £(S) 


il 
I 

ft  k 


400  M(1)». 0656/ .075 

420  M(2)*.303/.075 

440  M(3)»1./.075 

460  PRINT, "M2 fM3,M4",M 
480  A.1170. 

500  &i  .0656* 15600. 

520  C.  .303*15600.  ♦  .0656*175. 

540  D. 15600  .+0 .0656+ .303*175. 

560  E»  175 .303 

580  PRINT,"A,B,C,D",A,B,CfDfE 
600  pRiirr,"iNPur,N" 

610  INPUT ,N 

620  DO  400  IK»1,N 

640  PRINT, "R  COEFFS" 

660  INPUT, R2,R3,R4,R5 
680  PRINT,"R2,R3,R4,R5",R2,R3,R4,R5 
700  K4  «  R2-E 
720  K3  .  (B-R5J/B 
740  K2  *  (C*(R5/B)-R4)/B 
760  K1  •  (D*(R5/B)-C»K2-R3)/B 
780  PRINT, "K1 ,K2 ,K3  *K4" ,K1 ,K2 ,K3 ,K4 
800* 

820  IH2X-7 
825  H(1,7)*K2 
830  H(1,6)  ■  K3  +  K2*E 
840  H( 1 ,5)  ■  K4  +  K3*E  ♦  K2*D 
850  H(1,4)  •  -K1*B 
860* 

880  H(2,7)  *  K3 
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11:54:02 


SENSE  (continued) 


890  H(2,6)  •  K3*E+K4 
900  H(2,5)  ■  -K1*B  -  K2*C 
910  H(2,4)  ■  -K2*B 
920* 

940  H(3#7)  »  K4 

950  H<3,6)  «  -(K1*B  ♦  K2*C  ♦  K3»D) 
960  H<3,5>  *  -K2*B  -  K3*C 
970  H(3«4)  •  -K3*B 
980* 

990* 

1000  IF2X1*3 

1010  F( 1,3)  •  -K2 
1020  F(1,2)  *  -(K3-1.) 

1030* 

1050  F(2,3)  *  -U3-1.) 

1060* 

1079  IF2X2*4 

1080  F(3,4).1. 

1090  IDELA-3 
1091* 

1093  DELA1(3)*1. 

1094  DELA1(2)*.303 

1095  DELA1(1)».0656 
1097* 

1098  DELA2(3)»1. 

1099  DELA2<2)*175. 

1100  DEL A2( 1) *15600 . 

1102* 

1104  IDELF.5 

1105  DELF(5)*1 

1106  DELF(4).R2 

1107  DELF(3)«R3 

1108  DELF(2)*R4 

1109  DELF(1)*R5 
1110* 

1120  DO  100  1*1,2 

1130  100  1GN(I)*IF2X1 

1131  IGN( 3)*IF2X2 
1135* 

1140  DO  110  1*1,3 

1150  DO  110  J*1,4 

1160  110  GN(Z,J)  *  M(I)*F(I,J) 

1170* 

1180  IGN(4)«1DELA 
1190  IGN(5)-1DELA 
1200  DO  115  1*1,3 
1210  GN(4,I)*DELA1(1) 

1220  115  GN(5,l)aDELA2(l) 

1270* 

1280  DO  120  1*1,3 
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SENSE  (continued) 

1290  1GN(5+I)*IH2X 
1300  120  CONTINUE 
1310* 

1320  DO  130  1*1,3 
1340  DO  130  J*4,IH2X 
1350  130  01(5*1, J)  •  M(I)*H(I,J) 
1360* 

1370  KN(9)»1 
1380  CN(9,1)*1 
1390* 

1400  IGN( 10)*1 
1410  GN(10,1)*1170. 

1420* 

1430* 

1440* 

1450  IGD(1)«IDELA 
1460  IGD(2)*ID£LA 
1480  00  140  J*1,3 
1490  140GD(1,J)*DELA1(J) 

1500  DO  150  J*1,3 
1510  150  GD(2,J)»DBLA2(J) 

1520* 

1530  IGD(3)*IDELF 
1540  DO  160  1*1,5 
1550  160  GD(3,I)*DELF(I) 

1560* 

1570* 

1580  IGNP* 10 

1590  IGDP.3 

1591  IGDP*4 

1592  GD(4,2)*1 

1593  IGD(4)*2 

1600  P2(6)."A" 

1601  P2(7)*"*" 

1602  P2(8)*"f" 

1615* 

1620  DO  180  1*1,40 
1630  CV(1,I)*PKI) 

1640  180  CV(2,I)»P2(I) 

1645  PRINT, "PI", (P1(I), 1*1, 10) 

1646  PRINT, "PI" , (PI (I) ,1*11,20) 

1647  PRINT,"P2",(P2(I),I*1,10) 

1660  WN*7 

1661  WI*U)1 
1670* 

1680  CALL  SEN SUB 
1685  400  CONTINUE 
1690* 

1700  END 
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SENSEM 


100***  MAIM  PROGRAM  TO  COMPUTE  THE  SEN6ITIVITY 
120***  OF  E  VITH  RESPECT  TO  M(I) 

140*** 

160*** 

180  COMMON  GN,GD,IGN,IGD,WN,VI,IGNP,  I  GDP ,  CV 

185  COMMON  ICNT 

200  IMPLICIT  DOUBLE  PRECISI0N(A-H,0-Z) 

220  DOUBLE  PRECISION  K(4),M 

240  INTEGER  IGN< 20) , IGD( 20) , VN 

260  DIMENSION  GN( 20 , 10 ) ,GD( 20 , 1 0) 

280  DIMENSION  F< 4 , 4 , 4 ) , IF< 4 , 4 ) 

300  DIMENSION  M(4) ,DELA1 ( 3 ) ,DELA2(3) ,DELF( 7) 

320  CHARACTER  CV(2,40) 

340  CHARACTER* 1  PI (40) ,P2(40) 

360  DATA  PI /"A","#"/ 

380  DATA  P2/"A"  ,  "A"  ,  "*" ,  *'#"/ 

400  LIBRARY  "SENSUB" 

420  PI ( 2 
425  ICNT-0 
440* 

460* 

480* 

500  M(2)«. 0656/. 075 

520  M(3)-.303/.075 

540  M(4)-l./.075 

560  PRINT,"Ml,M2,M3,M4",M 
580  A-1170. 

1980  IDELA-3 
2000  DELA1  (  3  )- 1  . 

2020  DELA1 ( 2 )- . 303 
2040  DELA1 ( 1  )•  . 0656 
2060* 

2080  DELA2 ( 3  )- 1  . 

2100  DELA2 ( 2 )• 1 75 . 

2120  DELA2 ( 1 )- 1 5 600 . 

2270* 

2280*  LOOP  FOR  THE  THREE  M  VALUES 
2290* 

2300  DO  400  11*2,4 
2310  DO  300  1-1,10 
2320  300  GN( 1  , 1  )-0 . 0 
2330  IGN( 1 )-I 1-1 
2340  GN ( 1 , 1 1  - 1 )*M( II)*A 
2580  IGD( 1 )*IDELA 
2600  DO  140  J-l  ,3 
2620  140CD( 1 , J)-DELAl(J) 

2625  IGD( 2 )* IDELA 

2627  DO  150  1-1,3 

2628  150  GD(2,I)-DELA2(1) 

2760  IGNP-1 
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SENSEM 


(continued ) 


2780  IGDP-2 

2860  DO  180  1-1,40 

2880  CV(1 ,I)-P1(1) 

2900  180  CV(2,1)-P2(I) 

2920  PR1MT,"P1,,,(P1(1),  1-1,10) 
2960  PRIIT,"P2M,(P2(1) ,1-1 , 10) 
3060  WN-6 
3080  HI-. 001 
3100* 

3120  CALL  SENSUB 
3140  400  CONTINUE 
3160* 

3180  END 


-1 


SENSEAO 


100***  MAIM  PROGRAM  TO  COMPUTE  THE  FREQUENCY 
120***  RESPONSE  OF  E(S) 

130***  VITH  RESPECT  TO  OBSERVER  CAIN  AND  THETA 
140*** 

iso*** 

180  COMMON  CN,CD,IGN,ICD,WN,WI,ICNP,IGDP,CV 

200  IMPLICIT  DOUBLE  PREC1SI0N(A~H,0-Z) 

220  DOUBLE  PRECISION  K 1 ,K2 ,K3  ,K4 ,M 

240  INTEGER  IGN(20) , ICD<20) ,VN 

260  DIMENSION  GN( 20 , 10 ) , CD( 2 0 , 1 0) 

2 6 5  DIMENSION  F( 3 , 4 ) , H( 3 , 7 ) 

266  DIMENSION  M(3) , DELA l ( 3)  , DELA2 ( 3 ) , DELF( 7) 

280  CHARACTER  CV(2,40) 

281  CHARACTER* 1  Pi (40) ,P2(40) 

285  DATA  P2/"A" 

287  DATA  PI /2*"A" , . "A" , "*" , "#"/ 

300  LIBRARY  nSENSUBH 
340* 

360* 

380* 

400  M( 1 )-. 0656/. 075 

420  M(  2)-. 303/. 075 

440  M(3)-l./.075 

460  PRINT, MM2 . M3, M4".M 
480  A-1170. 

500  B-. 0656*15600. 

520  C-. 303*15600.  ♦  .0656*175. 

540  D- 15600. ♦O. 0656*. 303*175. 

560  E-175.303 

580  PRINT,"A.B.C,D",A,B,C,D 
600  PRINT, "INPUT, N" 

610  INPUT, N 

620  DO  400  IK-1  ,N 

640  PRINT, "R  COEFFS" 

660  INPUT, R2,R3,R4,RS 
680  PRINT ,"R2,R3,R4(R5"aR2,R3,R4,R5 
700  K4  -  R2-E 
720  K3  -  (B-R5)/B 
740  K2  -  (C*(R5/B)-R4)/B 
760  K 1  «  (D*(R5/B)-C*K2-R3)/B 
780  PRINT , "K 1 ,K2 , K3 , K4" ,K1 ,K2 ,K3,K4 
800* 

1090  IDELA-3 
1091* 

1093  DELA1 ( 3)-l . 

1094  DELA1 ( 2 )» . 303 

1095  DELA1 ( 1 )• . 0656 
1097* 

1098  DELA2C 3 )* 1 . 

1099  DELA2( 2 )•! 75 . 
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SENSEAO  (continued) 


1100  DELA2 ( 1 )■ 1 5600 . 

1102* 

1104  IDELP-5 

1105  DELF(5)-1 

1106  DELP(4)-R2 

1107  DELF(3)-R3 

1108  DELF(2)«R4 

1109  DELF( 1  )-R5 
1110* 

1130  IGN(l)-2 

1140  GN(1 ,1)-H(1)*K3 

1150  GM( 1 ,2)«M( 1)*K2+M(2)*K3 

1  390* 

1400  I CM < 3)-l 
1410  GR( 3 , 1 )-- 1170. 

1420* 

1430* 

1440* 

1456  IGN( 2 )*IDELA 
1480  DO  140  J-1,3 
1485  GM  <  2 , J )-- M( 3 )*DELA1( J) 

1490  140  CONTINUE 
1520* 

1530  IGD(1)-1DELF 
1540  DO  160  1-1 ,5 
1550  160  GD< 1 , I)-DELF(I) 

1560* 

1570* 

1580  IGNP-3 
1590  IGDP-I 
1615* 

1620  DO  180  1-1,40 
1630  CV(1,I)-P1(I) 

1640  180  CV(2,I)-P2(I) 

1645  PRINT, "PI”, ( F 1 < I ) ,1-1 ,10) 

1646  PRINT,,,P1,,,(P1(I)  ,1-11 ,20) 

1647  PRINT, ,,F2”,(P2(I), 1-1,10) 

1660  WN-7 

1661  WI-.01 
1670* 

1680  CALL  SEN5UB 
1685  400  CONTINUE 
1690* 

1700  END 
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M  RS"  ‘  PB'NTfD  IN  U  5  A 


-1- 

w  SENSEAS 

100***  MAIN  PROGRAM  TO  COMPUTE  THE  FREQUENCY 
120***  RESPONSE  OF  E(S) 

130***  WITH  RESPECT  TO  SENSOR  GAIN  AND  THETA  DOT 

140*** 

160*** 

180  COfMON  GN,GD,IGN,IGD,WN,NI,IGNP,IGDP,CV 

200  IMPLICIT  DOUBLE  PRECISION(A-HfO*Z) 

220  DOUBLE  PRECISION  K1 ,K2,K3,KN,M 

240  INTEGER  IGN(20),IGD<20),WN 

260  DIMENSION  GN(20,10),GD(20,10) 

265  DIMENSION  F(3,4),H(3.7) 

266  DIMENSION  M(3),DELA1(3),DELA2(3>,DELF<7) 

280  CHARACTER  CV(2,40) 

281  CHARACTER* 1  PI (40) ,P2(40) 

285  DATA  P2/2*"A" tnA" t"6*/ 

287  DATA  P1/2*" A" , , 2*" A" ,  , "A" , , "A” , *** , "€"/ 

300  LIBRARY  "SENSUB" 

340* 

360* 

380* 

400  M(1)«  .0656/ .075 

420  M(2)b .303/ .075 

440  M(3)«1./.075 

460  PRINTf"M2fM3,M4,,lM 
480  A»1170. 

500  Bb  .0656* 15600 . 

520  C. .303*15600.  ♦  .0656*175. 

540  D.  15600  ^0 .0656*  .3Q3M75 . 

560  Eb  175.303 

580  PUNT,"A,B,C,D*',A,B,CID,E 
600  PRUrr.-INPUT.N" 

610  INPUT, N 
620  DO  400  IKbI.N 
640  PRINT, "R  COEFFS" 

660  INPUT, R2,R3,R4,R5 
680  PRINT,NR2,R3,R4,R5N,R2,R3»R4,R5 
700  K4  •  R2-E 
720  K3  •  (B-R5)/B 
740  K2  •  (C*(R5/B)-R4)/B 
760  K1  -  (D*(R5/B)-C*IC2-R3)/B 
780  PRINT, NK1,K2,K3 ,K4",K1 ,K2,K3,K4 
800* 

1090  IDELA*3 
1091* 

1093  DCLA1(3)*1 . 

1094  D£LA1(2)»  .303 

1095  DELA1(1)b.0656 
1097* 

1098  DELA2( 3)b1 . 

1099  DELA2(2)b175. 
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^  SENSEAS  (continued) 

1100  DELA2(  0*15600. 

1102* 

1104  1DELF«5 

1105  DELF(5)a1 

1106  DELF(4)*R2 

1107  DELF(3)*R3 

1108  DELF(2)aR4 

1109  DELF(1)*R5 
1110* 

1130  IGN(1)*2 

1140  GN< 1,1 )*M(1)«K3 

1150  GN( 1 ,2)*M( 1 )*K2+M(2)*K3 

1390* 

1400  IGM<6)*1 
1410  GN(6, 0*1170. 

1420* 

w  1430* 

1440* 

1450  IGD(1)*1 

w  1455  IGN(3)*IDELA 

1456  IGN(2)*IDELA 

1457  IGN(4).IDELA 

W  1460  IGD(2)aIDELA 

1480  DO  140  J»1,3 
1485  GN(3,J)*-H(3)*DELA1(J) 
w  1490  140  CONTINUE 

1495  GD( 1 , 1 )■ 1 . 

1500  DO  150  J.1,3 
1505  GN(2,J)*DELA2( J) 

1507  GN(4,J)»DELA2(J) 

1510  150  GD(2,J)*DELA2(J) 

1520* 

1530  IGD(3)-IDELF 
1535  IGN(5)*IDELF 
*!  1540  DO  160  1*1,5 

1545  GN(5,I)»M(3)*DELF(I) 

1550  160  GD(3,I)*DELF(I) 

4  1560* 

1570* 

1580  IGNP*6 
'T '  1590  IGDP.3 

!  1615* 

1620  DO  180  1*1 ,40 
V  1630  CV(1,I)*P1(I) 

1640  180  CV(2,I)*P2(I) 

1645  PRINT, "PI", <P1(I), 1*1, 10) 

^  1646  PRINT, "PI",  (P1U)  ,1*11,20) 

■r  1647  PRMT,"P2",(P2(I) ,1*1,10) 

1660  HN*7 

^  1661  HI*  .01 

124 


07  FEB  79 


vs  r.  Nt  a-UfJKJd  oi ss  vv 


SEN SEAS  (continued) 


*1* 


1670* 

1680  CALL  SENSUB 
1685  NOO  CONTINUE 
1690* 

1700  END 


<- 

r  r*-* 


s' 
*  _ 

T. 


A* 
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1  FORM  8610  Pk  MIkD  IN  USA, 


1 


SENSEZ1 

100***  MAIN  PROGRAM  TO  COMPUTE  THE  FREQUENCY 
120***  RESPONSE  W/R  TO  ZETA  1 

140»*» 

160*** 

180  COhMON  GN,GD,IGN,IGD,WN,WI,IGNP,IGDP,CV 

181  OOIMON  ICNT 

200  IMPLICIT  DOUBLE  PRECISION (A*H,0*Z) 

220  DOUBLE  PRECISION  K1 ,K2,K3,K4,M 

240  INTEGER  IGN(20) ,IGD(20) ,WN 

260  DIMENSION  GN(20,10) ,GD(20,10) 

280  DIMENSION  CC(4,7) 

300  DIMENSION  M<4) ,DELA1(3) ,DELA2(3),0ELF(7) 

320  CHARACTER  CV<2,40) 

340  CHARACTER* 1  P1(4Q) ,P2(40) 

360  DATA  P2/3•"A%2•,,•^2«nV*\"*,,, "A", 

380  DATA  P1/3*"A",2*%*,"A","*,‘,2*"A,'t2*"*",’'§"/ 
400  LIBRARY  "SENSUB" 

402  ICNT*0 
420  P2(4)»"*" 

440  P2(5)«"*" 

460* 

480* 

500* 

520  M(2)«. 0656/ .075 

540  M(  3)*. 30  V  .075 

560  M(4).1./.075 

580  PRINT, "Ml  ,M2,M3,M4,,,M 
600  A.1170. 

620  B.  .0656*  15600. 

640  C. .303*15600.  +  .0656*175. 

660  D*  15600  .♦0.0656+ .303*  175. 

680  E*  175 .303 

700  PRINT,"A,B,C,D,,,A,B,C,D,E 
705  PRINT, "INPUT, N" 

710  INPUT, N 
720  DO  400  1.1, N 
722  PRINTER  COEFFS" 

740  INPUT, R2,R3,R4,R5 
800  PRINT,"R2,R3*R4,R5,,*R2rR3»R11»R5 
820  K4  .  R2-E 
V '  840  K3  *  (B-R5)/B 

860  K2  «  (C*(R5/B)-R4)/B 
880  K1  .  (D»(R5/B)-C*K2-R3)/B 
v'  900  PRINT, HK1  ,K2,K3,K4",K1  ,K2,K3»K4 

920* 

940  IH2X.7 

w  960  CC(2,7)-K2 

980  CC(2,6)  •  K3  ♦  K2*E 
1000  CC(2,5)  ■  K4  ♦  K3*E  ♦  K2*D 
w  1020  CC(2,4)  ■  -K1*B 
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1040* 

1060  CC(3,7)  »  K3 
1080  CC(3,6)  •  K3*E+K4 
1100  CC(3*5)  •  -K1*B  -  K2*C 
1120  CC(3,4)  -  -K2*B 
1140# 

1160  CC(4,7)  ■  K4 

1180  CC(4,6)  .  -(K1*B  ♦  K2»C  ♦  K3*D) 
1200  CC(4,5>  *  -K2*B  -  K3*C 
1220  CC(4,4)  *  -K3*B 
1240* 

1260* 

1280  1DELA*3 
1300* 

1320  DELA1(3)»1. 

1340  DELA1<2)«.303 
1360  DELA1(1)*.0656 
1380* 

1400  DELA2(3>»1. 

1420  DELA2(2)*175. 

1440  DELA2(1)» 15600. 

1460* 

1480  IDELF*5 
1500  DELF(5)*1 
1520  DELF(4).R2 
1540  DELF(3)*R3 
1560  DELF(2)*R4 
1580  DELftl)»R5 
1600* 

1620* 

1640  IGN(1)*7 

1660  IGN(2)»7 

1680  IGN(3)»7 

1700  DO  100  1*2,4 

1720  DO  100  J*4,7 

1740  100  GN(I-1,J)*M<I)*CC<I,J) 

1760* 

1780  IGN(4)»1 
1800  GH(4,1)..5915*A 
1820* 

1840  IGN(5)*2 

1860  GN(5,2)*2.*SQRT(  .0656) 

1880* 

1900  IGN(6)*3 

1920  DO  110  1*1,3 

1940  110  GN(6,I)*DELA2(I) 

I960* 

1980* 

2000* 

2020* 


-> 


SENSEZ1  (continued) 

2040  IGD(  1)*ID£LA 
2060  IGD(2)*IDELA 
2080  IGD(3)*IDELA 
2100  IGD(4)aIDELA 
2120  DO  140  J*1,3 
2140  GD(2,J)*DELA1(J) 

2160  140GDC 1 ,J)*DELA1( J) 

2180  DO  150  J-1,3 
2200  GD(3»J)«DELA2(J) 

2220  150  GD(4,J)aD£LA2( J) 

2240* 

2260  ZGD(5)a2DELF 
2280  DO  160  1*1 ,5 
2300  160  GD(5,I)*D£LF(I) 

2320* 

2340  IGD(6)a2 
2360  GD(6,2)a1 
2380» 

2400  IGNP*6 
2420  IGDPa6 
2440* 

2460  DO  180  1*1,40 
2480  CV(1,I)«P1(I) 

2500  180  CV(2,I)»P2(I) 

2520  PRINT, "PI", (P1(I), 1*1, 17) 
2540  PRINT,"P2",(P2(I),I«1,17> 
2560  WN-6 
2580  Vil*.001 
2600* 

2640  CALL  SENSUB 
2645  400  CONTINUE 
2660  END 


* 
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1620* 

1640  IGN ( 1  )-  7 

1660  IGN (  2  )-7 

1680  IGN (  3  )"7 

1700  DO  100  1-2,4 

1720  DO  100  J-4,7 

1740  100  GNCl-1 , J)-M(l)*CC(ItJ) 

1760* 

1780  IGN(4)-1 

1800  GN(4,1)-A*SQRT( .0656) 

1820* 

1840  IGN( 5 )-2 

1860  GN(5,2)-2.*( .5915) 

1880  GM<5,1)-2.*SQRT( .0656) 
1900* 

1920  ICN( 6 )-3 

1940  DO  110  1-1,3 

1960  110  GN(  6 , I )-DELA2 ( I ) 

1980* 

2000* 

2020* 

2040* 

2060  IGD( 1 )-IDELA 
2080  IGD( 2 )>IDELA 
2100  IGD( 3 )-IDELA 
2120  IGD(4)-IDELA 
2140  DO  140  J-1,3 
2160  GD( 2 , J )-DELAl( J ) 

2180  140GD(1 ,J)-DELA1(J) 

2200  DO  150  J-1,3 
2220  GD( 3 ,  J )-DELA2(j ) 

2240  150  GD(4 , J ) -DELA2 ( J) 

2260* 

2280  IGD(5 )-IDELF 
2300  DO  160  1-1,5 
2320  160  GD(5, I)-DELF(I) 

2340* 

2360  IGD(6)*2 
2380  GD(6,2)-1 
2400* 

2420  IGNP-6 
2440  IGDP-6 
2460* 
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1620* 

^  1640  IGM(1).7 

1660  IGM(2)«7 
1680  1GH(3)*7 
1700  DO  100  1.2,4 
1720  DO  100  J.4,7 
1740  100  GN(I-1,J).M(I)*CC<I,J) 
1760* 

1780  IGN(4)»1 

1800  GN(4,1)«A»SQRT(  15600.) 

'  1820* 

1840  IGN(5)*2 
1860  GN(5,2)«2.* .70056 
^  1880  GM(5,1).2.*124.9 

1900* 

1920  IGN(6)«3 

^  1940  DO  110  1.1,3 

I960  110  GN(6,1).DELA1(I) 

1980* 

2000* 

2020* 

i.  2040* 

^  2060  IGWD.IDELA 

2080  IGD(2)«IDELA 
2100  IGD( 3).IDELA 
2^  2120  KD(4).IDELA 

!?  2140  DO  140  J.1,3 

<  2160  G0(2,J).DELA1(J) 

rV  2180  140GD(1,J).DELA1(J) 

2200  DO  150  J.1,3 
2220  GD(3,J)*DELA2(J) 

2240  150  GD(4,J).DELA2(J) 

2260* 

2280  IGD(5).1DELF 
i."  2300  DO  160  I«1,5 

2320  160  GD(5,I).DELF(I) 

2340* 

2360  IGD(6)>2 
2380  GD(6,2).1 
2400* 

"  2420  1GNP.6 

2440  IGDP.6 
2460* 


130 


07  FEB  79 


12:15:47 


SENSEZ2 


1640* 

1660* 

1680  IGN<1)-7 
1700  IGN<2)-7 
1720  IGN(3)«7 
1740  DO  100  1-2,4 
1760  DO  100  J-4,7 
1780  100  GNU-1, J)-M<I)»CC(I,J> 
1800* 

1820  IGN(4)«1 

1840  GN(4,1)-A* .700561 

I860* 

1880  IGN(5)-2 
1900  GN(5,2)-2.*124.9 
1920* 

1940  IGN(6)*3 

I960  DO  110  1-1,3 

1980  110  GN(6,I)«DELA1(I) 

2000* 

2020* 

2040* 

2060* 

2080  IGD(1)«IDELA 
2100  IGD(2)*IDELA 
2120  IGD(3)-1DELA 
2140  1GD(4)-IDELA 
2160  DO  140  J-1,3 
2180  GD(2VJ)-DELA1(J) 

2200  140GD(1tJ)»DELA1(J) 

2220  DO  150  J-1,3 
2240  GD(3,J)-DELA2(J) 

2260  150  GD(4,J).DELA2(J) 

2280* 

2300  IGD(5)*1DELF 
2320  DO  160  1-1,5 
2340  160  GD(5,I)*DELF(Z) 

2360* 

2380  IGD(6)-2 
2400  GD(6,2)-1 
2420* 

2440  IGNfr6 
2460  IGDP-6 
2480* 
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2260*** 

w  2262**  HERE  FOR  THE  VARIOUS  K  VALUES,  Hal,* 

2283* 

2300  Ilal 

W  2320*  00  *00  Ilal,* 

23*0  IGN(1)aIF(2,II) 

2360  DO  70  Ia1,IGN(1) 

V  2380  70  GN(1,I)aM(2)*F(2,II,I) 

2*00  IGN(2)-IF(3,II) 

2*20  00  80  Ia1,IGN<2) 

^  2**0  80  GN(2,I)aM(3)*F(3,H,I) 

2*60  IGII<3>aIF<*,II) 

2*80  DO  90  Ia1,IGN(3) 

"  2500  90  OI(3,I)aM(*)*F(*,II,I) 

2520  IGN(*)a* 

25*0  GN(*,*)aK(II)*A 
2560* 

2580  XGD(1)aIDELA 
2600  DO  140  Ja1,3 
2620  1*0GD(1,J)aDELA1(J) 

2625  IGD(4)aIDELA 

2627  DO  150  la  1,3 

2628  150  GD(*,I)»DELA2(I) 

26*0* 

2660  IGD(2)aIDELF 
2680  DO  160  1*1 ,5 
2700  160  GD(2,I)a0ELF(I) 

2720* 

27*0* 

2760  IGNPa* 

2780  IGDPa* 

2800  GD(3,2)a1 
2820  IGD(3)-2 
28*0* 
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100***  MAIN  PROGRAM  TO  COMPUTE  THE  SENSITIVITY 
120***  OF  ECAP  W/R  TO  ZETA1  OF  THE  OBSERVER 

140*** 

160*** 

180  COmON  GN,GD,IGN,IGD,WN,WI,IGNP,IGDP,CV 

185  COtMON  ICNT 

200  IMPLICIT  DOUBLE  PRECISION (A-H.O-Z) 

220  DOUBLE  PRECISION  K(4),M 

240  INTEGER  IGN(20)  .IGDC20)  ,WN 

260  DIMENSION  GN( 20, 10) ,00(20,10) 

280  DIMENSION  F<4,4,4) ,IF(4,4) 

300  DIMENSION  M(4),DELA1(3)>DELA2(3).DELF(7) 

320  CHARACTER  CV(2,40) 

340  CHARACTER* 1  P1<40) ,P2<40) 

360  DATA  P1/3*"A",2*"+",2*"A",2*"*","#"/ 

380  DATA  P2/1HA,1HA,1HA,1H*,1H*,1HA,1H*,1Hf/ 

400  LIBRARY  "SENSUB" 

425  ICNT*0 
440* 

460* 

480* 

500  M( 2) *.0656/ .075 

520  M( 3) *.30 3/ .075 

540  M(4)*1./.075 

560  PRINT, "Ml ,M2,M3,M4",M 
580  A*1170. 

600  B.  .0656*15600. 

620  C» .303* 15600.  ♦  .0656*175. 

640  D*  15600  .♦0.0656«- .303*  175. 

660  £*175.303 

680  PRINT,"A,B,C,D",AyB,C,D,E 
700  PRINT, "INPUT, N" 

710  INPUT,N 

720  DO  400  JJ*1,N 

730  PRINT, "R  COEFFS" 

740  INPUT, R2,R3,R4,R5 

780  PRINT, "R2,R3,R4,R5",R2,R3,R4,R5 

800  K(4)  ■  R2-E 

820  K(3)  «  (B-R5)/B 

840  K(2)  ■  (C»(R5/B)*R4)/B 

860  K(1)  *  (D*(R5/B)-C*K(2)*R3)/B 

880  PRINT, "K1 ,K2,K3,K4",K 

900  DO  10  1*1,4 

920  DO  10  J*1,4 

940  DO  10  KL*1,4 

960  10  F(I,J,KL)*0. 

980** 

1000  IF(2,1)*2 

1020  F(2,1,1)«B*(K(3)-1) 

1040  F(2,1,2)*B*K(2) 
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1060* 

1080  IF<2,2)«4 

1100  F(2,2,2)»-B»K(1)-0»(K(3)-1) 
1120  F(2,2,3)»E*K(4) 

1140  F(2,2»4)«1 
1160* 

1180  IF(2,3)*3 
1200  F(2f3,1)—K(1)*B 
1220  F(2,3,2)-K(2)«D  ♦  B*K<4) 
1240  F(2,3»3)»1 
1260* 

1280  IF(2f4)«3 
1300  F(2,4,2)«-K(3)+1 
1320  F(2,4,3)»-K(2) 

1340* 

1360  IF(3,1)*2 

1380  F(3,1,2)«B»(K(3)-D 

1400* 

1420  IF<3,2)*2 

1440  F(3,2f1)*B»(K(3)-1) 

1460  F(3»2,2)«C*(K(3)-1) 

1480* 

1500  IF(3,3)»4 
1520  F(3,3,D— B»K(2) 

1540  F(3,3,2).B*K(1)-C*K<2) 

1560  F(3,3,3)-K(4)>E 
1580  F(3,3,4)»1 
1600* 

1620  IF(3*4)*3 
1640  F<3,4,3)*1-K<3> 

1660* 

1680  IF(4,1)»3 
1700  F(4,1,3)— B 
1720* 

1740  IF(4,2)*3 
1760  F(4,2t2)»-B 
1780  F(4,2,3)*-C 
1800* 

1820  IF(4,3)«3 
1840  F(4,3t1).-B 
1860  F(4,3,2)—C 
1880  F(4,3.3)«-D 
1900* 

1920  IF(4,4)»4 
1940  F(4f4,4)«1 
I960* 

1980  I0ELA-3 
2000  DELA1(3)*1 • 

2020  DELA1(2)a.303 
2040  DGLA1(1)».0656 
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2060* 

2080  DELA2(3)»1. 

2100  DSLA2(2)»175. 

2120  DELA2(1)* 15600. 

2140* 

2160  10ELF>5 
2180  0ELF(5)>1 
2200  DELF(4)*R2 
2220  DELF(3)«R3 
2240  DELF(2)*R4 
2260  DELF(1)*R5 
2279*** 

2280*** 

2283* 

2300  11*4 

2340  IGN(1)*IF(2,II) 

2360  DO  70  I*1,IGN(1) 

2380  70  GN( 1 ,I)*M(2)*F(2,II,I) 

2400  IGN(2)*IF(3,II) 

2420  DO  80  I*1,IGN(2) 

2440  80  GN(2,I)»M(3)*F(3,II,I) 

2460  IGN(3)»IF(4,II) 

2480  DO  90  1*1 fIGN(3) 

2500  90  GH(3,D»M(4)*F(4,II.I) 

2520  IGN(4)*2 
2540  GH(4t2)*2*SQRT(.0653)*A#.5915 
2560* 

2580  IGD( 1)*1DELA 
2600  DO  140  J*1,3 
2620  140GD(1fJ)*DELA1(J) 

2625  IGD(4)*1DELA 

2626  1GN(5)*IDELA 

2627  DO  150  1*1,3 

2628  GD(4,I)*DELA2(I) 

2629  150  GN(5,I)*DELA2(I) 

2640* 

2660  IGD(2)*IDELF 
2680  DO  160  1*1,5 
2700  160  G1X2,I)*DELF(Z) 

2720* 

2740* 

2760  1GNP-5 
2780  IGDP*4 
2800  GD(3,2>*1 
2820  IGD<3)«2 
2840* 

2860  DO  180  1*1,40 
2880  CV(1,I)*P1(I) 

2900  180  CV(2,I)*P2(I) 

2920  PRWT,"P1*,(P1(I),I*1,10) 
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2960  PRMTf"P2",(I>2(I),I»1,10) 
3060  UN-6 
3080  W..001 
3100* 

3120  CALL  SEN SUB 
3140  400  CONTINUE 
3160* 

3180  END 
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2280*** 

2283* 

2300  1 1-4 

2340  IGM( 1 )-IF( 2,11) 

2360  DO  70  1-1 ,IG*(1) 

2380  70  GH<1 , I )-M( 2 )*F( 2 , I I , I ) 
2400  ZGR(2)-IF(3,II) 

2420  DO  80  1-1 ,IGM(2) 

2440  80  GH(2,I)-M(3)*F<3,II,I) 
2460  IGN(3)-IF<4,II) 

2480  DO  90  1-1 , IGM( 3) 

2500  90  GN(3,I)-M(4)*F(4,II,I) 
2520  lGN(4)-2 

2540  GH(4,2)-A*.303 

2541  GN(4,1)-2*A*.0656 
2560* 

2580  ICD(1)-IDELA 
2600  DO  140  J-1,3 
2620  140GD( 1 , J)aDELAl ( J ) 

2625  IGD(4)*IDELA 

2626  ICR( 5 )«IDELA 

2627  DO  150  1-1,3 

2628  GD(4 , I )-DELA2 ( I ) 

2629  150  GN(5,I)-DELA2(I) 

2640* 

2660  IGD( 2 )-IDELF 
2680  DO  160  1-1,5 
2700  160  GD( 2 , 1 )-DELF(  I ) 

2720* 

2740* 

2760  IGHP-5 
2780  IGDP-4 
2800  GD( 3 , 2 )- 1 
2820  IGD( 3 )*2 
2840* 
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2280*** 

2283* 

2300  II»4 

23«0  IGN( 1).IF(2,II) 

2360  DO  70  I«1,1GN(1) 

2380  70  GN(1,I)*M(2)*F<2,II,I) 
2400  IGN(2).IF(3,II) 

2420  DO  80  I«1,IGN<2) 

2440  80  GN(2,I)>N(3)*F(3«1Z»I) 
2460  IGN(3)*IF(4,II) 

2480  DO  90  I.1,IGN<3) 

2500  90  GN<3,I)»M<4)*F(4,II,I) 
2520  IGN(4)»2 

2540  GN(4,2).A*175 

2541  GN(4,1).2*A*15600 
2560* 

2580  IGD<1).IDELA 
2600  DO  140  J.1,3 
2610  GN(5,J)«DEU1(J) 

2620  140GD( 1 ,J)»DELA1( J) 

2625  IGD(4).IDELA 

2626  IGN( 5) .IDELA 

2627  DO  150  1.1,3 

2628  GD(4,I)«DELA2(I) 

2629  150  CONTINUE 
2640* 

2660  IGD(2)»I0ELF 
2680  DO  160  1.1,5 
2700  160  GD(2,I).DELF(I) 

2720* 

2740* 

2760  IGNP.5 
2780  IGDP.4 
2800  GD(3,2)>1 
2820  IGD<3)*2 
2840* 
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2280*** 

2283* 

2300  I 1*4 

2340  IGM(1 )*IF(2,II) 

2360  DO  70  1-1  ,  IGN ( 1  ) 

2380  70  CM(1 ,I)«M(2)*F(2,II,I) 

2400  IGM(2)-IF(3,II) 

2420  DO  80  1-1 , IGH( 2 ) 

2440  80  G1<2,I)-M<3)*F<3,II,I> 

2460  ICM(3)-IF(4,il) 

2480  DO  90  I-1,ZCN(3) 

2500  90  GM(  3 , I )-M(4 )*P( 4 ,  1 1 , 1 ) 

2520  IGll(4)-2 

2540  CM (4 , 2 )-2*A* . 700561*SQRT( 15600.) 
2560* 

2580  IGD(1)-IDELA 
2600  DO  140  J-1,3 
2610  GM(5,J)-DELA1(J) 

2620  140GD( 1 , J)«DELA1 ( J) 

2625  IGD(4 )“ IDELA 

2626  ICM ( 5 )*IDELA 

2627  DO  150  1-1,3 

2628  GD(4 , I)-DELA2( 1) 

2629  150  COMTIMOE 
2640* 

2660  IGD( 2 )-IDELF 
2680  DO  160  1-1,5 
2700  160  GD( 2  ,I)-DELF(I) 

2720* 

2740* 

2760  IGMP-5 
2780  IGDP-4 
2800  GD( 3 , 2 )- 1 
2820  IGD( 3 )  —  2 
2840* 
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100  SUBROUTINE  SENSUB 

120  OOHNON  GNUM,GDEM,IGN,IGD,WN,WI 

140  COMM C»l  IGNP.IGOP.C 

160  COttlON  ICNT,I0UT1 ,I0UT2 

180  IMPLICIT  DOUBLE  PRECISION A-H.O-Z) 

200  COMPLEX  P(20) ,S,XM1 ,XM2,R1 , R2,R3 

220  INTEGER  0UT,WN,IGN(20),IGD<20) 

240  REAL  PA,PB,MAG(200,2),PHAS(200,2) 

260  DIMENSION  GNUM(20,10),GDEM(20,10),WX(20),WS(21) 

280  DIMENSION  A1(20),A2(20),B1(20),B2(20) 

300  DATA  WS/1, 1.25, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,2, 2 .25, 2 .5,2.75,3,3.5,4, 

320  A  4. 5, 5, 6, 7, 8 ,9/ 

340  CHARACTER  YES,F1LE1 ,FILE2 

360  CHARACTER  FILE3,C(2,40) 

380  NCNT»21 

420  NSTEP-2 
600* 

620*  OUTPUT  TO  FILE  OR  TTY 
640» 

650  OUTaO 

660  IF(ICNT  .NE.  0)GO  TO  10 

680  5  CONTINUE 

700  PRINT, "OUTPUT  TO  A  FILE?" 

720  INPUT, YES 

740  IF  (YES  JIE.  "YES")  GO  TO  10 

760  PRINT, "FILE  NAME  IS  ?" 

780  INPUT, FILE2 

800  PRINT, "FILE  NAME  IS  ?" 

820  INPUT  ,FILE3 

840  OPENFILE  4,FILE3 

860  REWIND  4 

880  ENDFILE  4 

900  I0UT2-4 

920  OPENFILE  3.FILE2 

940  REWIND  3 

960  ENDFILE  3 

980  I0UT1«3 

1000  10  CONTINUE 

1020  ICNTal 

1040* 

1045  ICNTal 

1060  WRITE(OUT,610) 

1080  WRITE(0UT,620) 

1100  WRITE(0UT,640) 

1120  DO  45  Jal.IGNP 

1140  WRITE(0UT,650)(GNUM(J,I),I-1,Ia1,IGN(J)) 

1160  45  OOTINUE 
1180  WRT.£(0UT,660) 

1200  DO  55  J» 1 ,IGDP 

1220  WRITE(0UT,650)(GDEM(J,I),I-1,Ia1,IGD(J)) 
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1240  55  CONTINUE 

1260  WRITE (OUT, 850) 

1280  WRITE(0UT,610) 

1300  WRITE  (OUT, 800) 

1320  KlaO 

1340  NCWT2-0 

1360  100  DO  250  KZ«1,WN 

1380  DO  200  K*1,NCNT,NSTEP 

1400  W*W5(K)*WI 

1420  W2«-W*W 

1440  K1.K1+1 

1460  CALL  PQSD(A1 ,B1 ,IGNP,W2,GNUM,IGN) 

1480  CALL  PQSD( A2 ,B2 , IGDP , W2 ,GDEM , IGD) 

1500  DO  105  I-1,IGNP 

1520  105  A1(I)  ■  A1(I)*W 

1540  DO  107  1*1,  IGDP 

1560  107  A2(I)  *  A2(I)*W 

1580  PHAS0C1 ,1)*W 

1600  MAG(K1 , 1)»W 

1620  ISW*1 

1640  DO  108  1*1 ,IGNP 

1660  108  P(I)  *  CMPLX(SNGL(B1(I)) ,SNGL(A1(I))) 

1680  110  S*CMPLX(0.,0.) 

1700  J*0 

1720  1.0 

1740  120  1*1+1 

1760  125  J.J+1 

1780  IF(C(ISWfJ)  .EQ.  "A")  GO  TO  130 

1800  IF(C(ISW,J)  .EQ.  "*")  GO  TO  135 

1820  IF(C(ISW,J)  .EQ.  "♦")  GO  TO  140 

1840  IF(C(ISW,J)  .EQ.  "-")  GO  TO  145 

1860  IF(C(ISW,J)  .EQ.  "S")  GO  TO  150 

1880  IF(C(ISW,J)  .EQ.  "R")  GO  TO  155 

1900  IF(C(ISW,J)  .EQ.  "t")  GO  TO  160 

1920  PRINT, "UNKNOWN  CHARACTER", C(ISW,J) 

1940  STOP 

I960* 

1980  130  R3*R2 

2000  R2.R1 

2020  R1aP(I) 

2040  GO  TO  120 

2060* 

2080  135  R1«R1*R2 

2100  137  R2.R3 

2120  GO  TO  125 

2140* 

2160  140  R1.R1+R2 

2180  GO  TO  137 

2200* 

2220  145  R1.R1-R2 
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*  SENSUB  (continued) 

2240  00  TO  137 

2260* 

2280  150  Ssfil 

2300  GO  TO  125 

2320* 

2340  155  R3-R2 

2360  R2»R1 

2360  R1«S 

2400  GO  TO  125 

2420* 

2440  160  IF(ISW  .EQ.  2)  GO  TO  170 

2460  ISW*2 

2480  DO  167  I»1,IGDP 

2500  167  P(I>  ■  CMPLX(SNGL(B2d)),SNGL(A2(I))) 

2520  XMlsRI 

2540  GO  TO  110 

^  2560» 

2580  170  XM2-R1 

2600  CALL  PHASE(  REAL(  XM1 ) ,  AIHAGC  XM1 ) , SMGL(W) ,  PA) 

^  2620  CALL  PHAS£(REAL(XM2),AlMAG(XM2),SNGL(Vf),PB) 

2640  PHASCK1 ,2)  *  (PA-PB)»180./3.141593 

2660*  XM1*CMPLX( REALC XM1 ) , AIMAG( XM1 )*SNGL( W) ) 

^  2680*  XH2»CMPLX(REAL(XH2)fAlMAG(XM2)*SNGL(W)) 

2700  MAG(K1,2)  ■  CABS(XM1)/CABS(XH2) 

2720  WRIT£(OUT ,810)W,MAG(K1 ,2) ,PHAS(K1 t2) 

^  2740  NCNT2«NCNT2+1 

2760  IF(NCNT2  XE.  44)  GO  TO  200 

2780  NCNT2*0 

r  ■'  2800*  WRITE(0UT,670) 

2820*  WRITE«XJT,610) 

2840*  WRIT£(0UT,800) 

2860  200  CONTINUE 
2880*  WRITE (OUT ,820) 

2900  WI-WI*10 

2920  250  CONTINUE 

2940  IFOICMT  XQ.  1)  GO  TO  270 

2960  NCNT»1 

2980  WN-1 

3000  GO  TO  100 

3020  270  CONTINUE 

1  3040  NCNT2*(55*NCNT2)/3 

3060  IF(I0UT1  .EQ.  0)GO  TO  320 

3080  DO  300  1.1, K1 

w  3100  300  WRITE( I0UT1 ,735)MAG( 1,1) ,HAG(  1,2) 

3120  DO  310  1.1, K1 

3140  310  WRITE( I0UT2 ,735) PHAS( 1,1), PHASC I ,2) 

w  3160  WRITE(I0UT1,745) 

3180  WRITE(I0UT2,745> 

3200  WRITECOUT ,850) 

w  3220  320  CONTINUE 
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3225  RETURN 
3240  • 

3260  * 

3280  • 

3300  • 

3320  600  FORMAK/,  3Q<"*"),  "INPUT  DATA",30("»")/) 

3340  690  FORMAK  V) 

3360  735  FORMAT(E13.5,"  ,  ",E12.5) 

3380  745  FORMAK"  1E37  1E37") 

3400  760  FORMAK 1H  ,"NU«ER  OF  DECADES  AND  INITIAL  FREQUENCY") 

3420  700  FORMAK 1H  ."GAIN  OF  G(S)*") 

3440  710  FORMAK  1H  ,*H0W  MANY  PXYNOMINALS  IN  THE  NUMERATOR  OF  G(S)") 

3460  720  FORMAK  1H  ,"HOM  MANY  PQLYNOMINALS  IN  THE  DENOMINATOR  OF  G(S)") 

3480  610  FORMAK 1H  ,70<"*")///) 

3500  620  FORMAK 1H  ,"G(S)") 

3520  640  FORMAK 1H  ."NUMERATOR  TERMS") 

3540  660  FORMAK 1H  ."DENOMINATOR  TERMS") 

3560  670  FORMAK/) 

3580  650  FORMAK 1H  ,4<E11 .5,"*S**",I2,"+")/5(5X,3(E11 .5,  "*S*«",  12,  "♦")/)) 
3600  800  FORMAK 1H  , "FREQUENCY", 3X,"MDB  OP-LOOP", 3X, "PHASE  0P-LP",3X, 

3620  4) 

3640  810  FORMAK  IX, G8 «3*G14.6,G14.6) 

3660  820  FORMAK  1H  ) 

3680  850  FORMAK////) 

3700  END 

3720* 

3740* 

3760* 

3780  SUBROUTINE  TEREAD(Y,J,IN,M,L1,C) 

3800  DOUBLE  PRECISION  Y(2Q,10) 

3820  DIMENSION  J(20) 

3840  CHARACTER  C(2,40) 

3860  READ(IN,999)M 

3880  IF(M  .EQ.  0)G0  TO  70 

3900  HRITE(OUT,998) 

3920  DO  40  1*1, M 

3940  READ(IN,999)Kf(Y(I,K+1-L),L*1,K) 

3960  J(I)  *  K 

3980  40  CONTINUE 

4000  IF(M  .EQ.  DGO  TO  60 

4020  WRITE(OUT ,997) 

4040  READ(IN,999)(C(L1,I) ,1*1,40) 

4060  45  OONKNUE 

4080  999  FORMAT(V) 

4100  997  FORMAK"  ENTER  THE  CHARACTER  STRING  -  ONE  LINE") 

4120  998  FORMAK  1H  ."ENTER  NWBER  OF  COEFFS  THEN  00EFFS,HIGH",/,5X, 

4140  &  "  TO  LON,  INCLUDING  ZEROS  FOR  EACH  POLYNONINAL”) 

4160  50  RETURN 

4180  60  C(L1,1)*"A" 

4200  C(L1,2)*"r 


SEN SUB  (continued) 


4220 

GO  TO  50 

4240  70 

Mai 

4260 

Y<1,1)»1. 

4280 

GO  TO  60 

4300 

END 

4320* 

4340* 

4360* 

4380 

SUBROUTINE  PHASE(R,I,W,P) 

4400 

REAL  I 

4420 

PI*3. 14 15926 

4440  10 

IF(  1)20,30 ,40 

4460  20 

IF(R)100,12Q,130 

4480  30 

IF( R) 100,140 ,140 

4500  40 

IF(R) 170,150,160 

4520  • 
4540  100 

PkPI+ATAN(I/R) 

4560  110 

RETURN 

4580  120 

P»PI«*.pI/2 . 

4600 

GO  TO  110 

4620  130 

P*2*PI+ATAN( I/R) 

4640 

GO  TO  110 

4660  140 

PmO. 

4680 

GO  TO  110 

4700  150 

PnPI/2. 

4720 

GO  TO  110 

4740  160 

P»ATAN(I/R) 

4760 

GO  TO  110 

4780  170 

P»PI+ATAN(I/R) 

4800 

GO  TO  110 

4820 

END 

4840* 

4860* 

4880* 

4900  SUBROUTINE  PQSD(A,B,M,Q,X,ID) 

4920  DOUBLE  PRECISION  X(20,10) ,A(20) ,B(20) ,Q 

4940  DIMENSION  ID(20) 

4960  DO  3  K-1.M 
4980  A(K)«0 
5000  B(K)*0 
5020  J*ID(K) 

5040  1  IF(J)  3,3,2 
5060  2  Z«B(K) 

5080  B(K)*Q»A(K)+X(M) 

5100  A(K)»Z 
5120  J.J-1 
5140  GO  TO  1 
5160  3  CONTINUE 
5180  RETURN 

5200  END 
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% 

POLYPAK 

1*  POLYPAK 

% 

2* 

3* 

4* 

POLYPAK  IS  A  PACKAGE  OP  ‘FORTRAN*  SUBROUTINES  FOR  POLYROMIRAI 

V 

5*  OPEKATIONS.  E.E.  MITCHELL  9/78 

6* 

7* 

ALL  COEFFICIENT  VECTORS  RAVE  A  DIMENSION  OF  20. 

s 

8* 

9* 

EACH  VECTOR  OF  COEFFICIENTS  MUST  BE  ORDERED  FROM 

***  LOW  ORDEI 

10*  TO 

HIGH  ORDER***  AMD  WITH  EACH  VECTOR  OF  COEFFICIENTS  THERE  MUSI 

11*  be 

AN  INTEGER  TELLING  THE  *NUMBER  OF  COEFFICIENTS* 

IN  THE  POLY. 

12* 

13* 

- 

14* 

15* 

THE  USAGE  IS  AS  FOLLOWS: 

16* 

'w 

17*  CALL  PADD(X,IX,Y,IY,Z,IZ) 

18* 

POLY  Z  IS  ADDED  TO  POLY  Y  AMD  GIVES  POLY  X 

19* 

20*  CALL  PSU>(Z,IZ,X,IX,Y,IY) 

21* 

POLY  Y  IS  SUBTRACTED  FROM  POLY  X  -  RESULTS  IN  POLY  Z 

22* 

V 

23* 

24*  CALL  PMPY(Z,IZ,X,IX,Y,IY) 

25* 

POLY  X  IS  MULTIPLIED  TIMES  POLY  Y  TO  GIVE  Z 

V 

26* 

100*  * 

************ 

*  * 

110*** 

THIS  PACKAGE  IS  A  COLLECTION  OF  FORTRAN  CALLABLE 

V 

120*** 

SUBROUTINES  FOR  POLYNOMINAL  OPERATIONS 

130* 

140* 

150 

SUBROUTINE  PADD(X , IX , Y , IY , Z , IZ ) 

160*  ADD 

TWO  POLYS 

170 

DIMENSION  X(20),Y(20),Z(20) 

180 

DO  10  1-1,19 

190  10 

X(I)-0. 

200 

DO  20  I-l.IY 

V 

210  20 

X( I )-Y( I ) 

i 

220 

DO  30  I-l.IZ 

230  30 

X(I)-X(I)*Z(I) 

V 

240 

IF(IY-IZ)40, 40,60 

250  40 

IX-IZ 

260  SO 

RETURN 

Nw 

270  60 

IX-IY 

280 

GO  TO  50 

290 

END 

<0> 

300*  * 

***** 

.  - 

310* 

320 

SUBROUTINE  PSUB< Z , IZ , X , IX , Y , IY) 

W  . 

3  30* 

2 


POLYPAK  (continued) 

340*  SUBTRACT  ORE  POLY  PROM  ANOTHER 


V 

350 

DIMENSION  X(20)  ,Y(20) ,Z(20) 

360 

N"  IX 

370 

IP(IX-IY)  10,20,20 

V 

380 

10  N-IY 

390 

20  IP(  N  )  90,90,30 

400 

30  DO  80  I  -l,N 

V 

410 

IP<  I-IX  )  40,40,60 

420 

40  IF(  I-IY  )  50,50,70 

430 

50  Z(I)  -  X(I)  -  Y(I) 

V 

440 

GO  TO  80 

450 

60  Z(I)  -  -Y(I) 

460 

GO  TO  80 

470 

70  Z(I)  -  X(I) 

480 

80  CONTINUE 

490 

90  IZ  -  N 

500 

RETURN 

510 

END 

520* 

530*  ****** 

S40  SUBROUTINE  PMPY(Z , IZ ,X , IX ,Y, IY) 
550* 

s"  560*  MULTIPLIES  TWO  POLYS 

570  DIMENSION  Z( 20) ,X( 20) ,Y(20) 

580  IZ-IX+IY-1 
w  590  DO  30  I-l.IZ 

600  30  Z(I)-0. 

610  DO  40  I-l.IX 
w  620  DO  40  J-l.IY 

630  X-I+J-l 

640  40  Z( K)aX( I )*Y( J)+Z( X) 
v  650  RETURN 

660  END 
670* 

V  680*  ****** 

690* 

'•T 

*  • 
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APPENDIX  C.  COMPUTER  PROGRAMS  FOR  THE  TIME 

RESPONSE  OF  THE  SENSITIVITY  EQUATIONS 

In  Appendix  B,  programs  were  listed  which  generated  the 
polynomials  of  the  sensitivity  equations  for  each  of  the 
parameters . 

Several  very  small  programs  were  then  written ,  a  typical 
one  being  STRAJP,  which  is  listed  below. 


STRAJP  09/24/78  19:20:58 


100***  MAIN  PROGRAM  TO  COMPUTE  THE  TIME  EQUATIONS 
110***  FOR  THE  SENSITIVITY  OF  E  WITH  RESPECT 

120***  TO  THE  PARAMETERS  Zl,  Z2,  Wl,  W2 

200* 

400  LIBRARY  "TRAJEQ" 

478  PRINT,  "****" 

479  PRINT,  "****" 

480  PRINT,  "RUN  FOR  PARAMETER  Zl**" 

481  PRINT,  "***" 

2630  C(1,40)-"Z1**" 

2635  C (1 , 39) ""DTRAJZl" 

2640  CALL  TRAJEQ 

These  small  programs  were  overlaid  on  the  sensitivity 
programs  of  Appendix  B.  The  combined  program  then  called 
subroutine  TRAJEQ  which  output  the  polynomial  coefficients 
into  a  file.  The  program  INTEGX  would  then  read  the 
coefficients  from  the  file  and  run  a  time  solution. 

The  subroutine  TRAJEQ  and  program  INTEGX  are  listed  next. 


1 


TRAJEQ 


100  SUBROUTINE  TRAJEQ 

120  COMMON  GMUM,CDBM,1GN,XGD,HN,HI 

140  COMMON  ICNP , XCDP , C 

160  COMMON  XCNT , X0UT1 , 10UT2 

180*  IMPLICIT  DOUBLE  PRECISION( A-H ,0-E) 

200  XNTECER  IC( 20) , IR3 , XR2 , XR1 , IS 

220  INTEGER  OUT ,HN , ICI( 20 ) , IGD( 20) 

240  REAL  PA,  PB ,  MAG( 200,2), PHAS( 200,2) 

260  DIMENSION  CNUM( 20 , 10 ) ,GDEM( 20 , 10) ,WX( 20) , WS( 2 1 ) 

280  DIMENSION  S(20),  Al ( 20) ,A2( 20) , Bl ( 20) , B2 ( 20 ) 

300  DIMENSION  CX< 20 , 1 0) , R1 ( 20 ) , R2( 20) ,R3( 20) 

320  LIBRART  "L . ES*** : POLYPAE" 

340  CHARACTER  TES .PILEl , PXLE2 

360  CHARACTER  FILE3 ,C< 2 ,40) 

380* 

400*  OUTPUT  TO  FILE  OR  TIT 
420* 

440  OUT-0 

460  IF(ICNT  .HE.  0)60  TO  10 

480*  OUT-0 

SOO*  3  CONTINUE 

320*  PRINT , "OUTPUT  TO  A  FILE?" 

340*  INPUT, TES 

360*  IF( YES  .ME.  "TES")  GO  TO  10 

380*  PRINT, "FILE  NAME  IS  ?" 

600*  INPUT, FILE2 

620  OPENF1LE  4  ,C( 1 , 39) , "NUMERIC" 

640  IOUT1  -  4 

660  REMIND  4 

680  ENDFILE  4 

700  10  CONTINUE 


440  OU 

460  IF 

480*  OU 

500*  5 

520*  PR 

540*  IN 

560*  IP 

580*  PR 

600*  IN 

620  OPE 

640  10 

660  RE 

680  EH 

700  10  CO 

720* 

740  ICNT-1 
760  HR 

780  HR 

800  HR 

820  DO 

840*  HR 

860  45  CO 

880  HR 

900  DO 

920*  HR 

940  55  CO 

960  HR 

980  HR 

1000  I8H-1 
1020  Xtl-1 
1040  XR2-1 
1060  X83-1 
1080  11 ( 1 )— 1 


HRITE(OUT ,610) 

URITE(OUT ,620) 

HRITE(OUT,640) 

DO  45  J- 1 , IGNP 

HR ITE ( OUT ,630)(GNUM( J, I), I- 1,1-1, ICN( J) ) 

CONTINUE 

NRITE(OUT ,660 ) 

DO  55  J-l.ICDP 

HRITE(OUT,650)(GDEM(J,X) ,1-1 ,1-1 , ICD(J)) 

CONTINUE 

HR ITE (OUT, 8 50) 

HR ITE (OUT, 610) 
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1100  B2(l)al 
1120  B3(l)al 
1U0  00  108  I-l.XCNP 
1160  IG( X )  -  ICll(I) 

1180  DO  108  Jal , IG( I ) 

1200  108  GX( 1 , J)  -  GMUN( I , J ) 


1220 

no 

CONTINUE 

1240 

J-0 

1260 

1-0 

1280 

120 

I-I+l 

1300 

125 

J-J+l 

1320 

IF( C( I8W , J ) 

.EQ.  WAM) 

GO 

TO 

130 

1340 

XF(C( ISV, J) 

.EQ.  "*•*) 

GO 

TO 

135 

1360 

IF(C( ISV, J) 

.EQ.  "♦") 

GO 

TO 

140 

1380 

IF(C( ISV, J) 

.EQ. 

GO 

TO 

145 

1400 

IF(C(ISW, J) 

.EQ.  "8" ) 

GO 

TO 

ISO 

1420 

IF( C( ISV , J ) 

.EQ.  "B") 

GO 

TO 

155 

1440 

IF(C(ISV,J) 

•EQ.  "§“> 

CO 

TO 

160 

1460 

PRINT , "UNKNOWN  CHARACTER1 

»,C(ISW,J> 

1480 

STOP 

1500* 

1520 

130 

IB3-IB2 

1540 

DO 

400  E-1,IK3 

1560 

400 

B3(K)  -B2(E) 

1580 

XB2 

-  IR1 

1600 

DO 

405  EM.XB2 

1620 

405 

B2(K)  *  B1(K) 

1640 

IB1 

-  IG( I ) 

1660 

DO 

410  Kal ( IB1 

1680 

410 

B1 ( K)  -CX(I,E) 

1700  GO  TO  120 
1720* 

1725* 

1740  135  CALL  PMPY( S , IS , B1 , IB1 , B2 , IB2 ) 
1760  DO  136  K*1 , 18 
1780  136  B1(K)-8(B) 

1800  IBl-IS 
1820  137  IB2  -  IB3 
1840  DO  420  K-1.XB2 
1860  420  B2(E)  -  B3(K) 

1880  GO  TO  125 
1900* 

1920* 

1940* 

1960  140  CALL  PADD(B1 ,111 ,B1 ,1*1 ,B2,I*2> 
1980  GOTO  137 
2000* 

2020* 

2040  14$  CALL  PSUB(B1 ,191 ,B1 ,IB1 ,B2,I«> 
2060  GO  TO  137 
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TKAJIQ  (continued) 


2080* 

2083* 

2100  130  IS-IR1 
2120  00  430  K*1 ,  IB 
2140  430  S(K)  -  81 (K) 

2160  00  TO  123 
2180* 

2200* 

2220  133  183*182 
2240  00  440  K- 1,183 
2260  440  83(X)*B2(X) 

2280  182  •  181 
2300  00  450  KM, 182 
2320  450  82(X)  -  81(K) 

2340  181  -  IS 
2360  DO  460  K-1,181 
2380  460  81 (X)  -  8(K) 

2400  60  TO  125 
2420* 

2440* 

2460* 

2480  160  IF(ISW  .EQ.  2)  GO  TO  170 
2500* 

2320  DO  167  IM.ICDF 
2540  1G( I )*1GD( 1 ) 

2560  DO  167  J-l.IG(l) 

2380  167  GX( 1 , J )  •  GDEM(I,J) 

2600  WRITE (OUT ,640) 

2620  WRITE (OUT, 645) 

2640  168  WRITE(0UT ,650) (81 ( 1),1~1, 1*1,181) 

2649  PPP-IRl*. 00001 

2650  WRITEdOUTl)  PPP, ( 81 ( I ) , IM , 181 ) 

2660  ISW-ISW+1 

2680  GO  TO  (110, 110, 172), ISW 
2700  170  WRITE (OUT, 660) 

2720  WRITE (OUT ,645 ) 

2740  GO  TO  168 

2760  172  WRITE (OUT ,610) 

2761  RBTURR 
2780  * 

2800  * 

2820  * 

2840  600  F0RMAT(/,30("*") ,MIUPUT  DATA" ,30("-") /) 

2860  690  FORMAT(V) 

2880  735  fORMAKBl 3. 5 , "  ,  H,E12.5) 

2  900  745  F0RNAT("1 .E38  ,  1.B38") 

2920  610  FORMAT ( 18  ,70<"-")///) 

2940  620  FORMATdB  ,"C(8)"> 

2960  640  F0RMAT(1H  , "NUMERATOR  TERMS") 


2980  645  FORMATC”  COMBINED") 
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TEAJEQ  ( 


d) 

3000  660  F0RMAT(1H  , "DENOMINATOR  TEENS") 

3020  670  F0BMAT( /) 

3040  650  TORNATdH  ,4(El  3 . 7  , "*8**" ,  12  ."♦">/5<  SX ,3(E1 3 . 7  ,"*8**" ,  12 
3060  820  rOEMAT(lB  ) 

3080  850  70ENAT(////) 

3100  END 


M 
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1NTEGX 


1*  IN  TEG 
2' 

3*  DESCRIPTION— Solves  a  act  of  N  first-order  differential  equations 
4*  of  the  fora; 

10* 

20*  dXU)  -  G(I)  for  I  a  1,2, . ,N, 

30* 

40*  when  the  initial  conditons  are  known,  by  the  fourth-order  Runge  Kutta 
SO*  technique  and  provides  a  plot  of  T  vs  X(N1),  X(N2),  X(N3),  and 
60*  X(N4)  in  that  order  depending  on  selection  in  DATA  statement . 

70* 

80*  INSTRUCTIONS— The  functions,  G(I),  are  to  be  inserted  after  line  1290 
90*  and  before  line  1320  in  the  following  fora: 

1295  G(1)  ■  WHATEVER  IT  IS 

1300  G(2)  •  WHATEVER  IT  IS 

The  entire  aatheaatical  aodel  should  go  in  this  section  including 
algebraic  constraints. 

1330  0ATH,T2,T1 ,M 
1350  DATA  H1,N1,N2,N3,N4 

WHERE  H  a  Increaant  size 

T2  a  Initial  value  of  T 
T1  a  Final  value  of  T 
M  a  NuWber  of  variables  to  be  plotted 
HI  a  No.  of  inoreaants  of  H  between  print  outs  (oust  be  an 
Integer.) 

N1,N2,N3,N4  Identify  the  dependent  variables  to  be  printed  out. 


T  a  Independent  variable 
G(I)  a  Known  functions  a  dX(I)/dT 

•  •  •  •  HAIR  PROGRAM  •  •  •  • 

361  DIM  Q(20),D(20) 

362  PRIIfT*INPUT  DATA  FILE  NAME" 

363  INPUT  F$ 

364  FILEil:  F$ 

365  READil:  17 

366  FOR  la  1  TO  17 

367  READil:  (XI) 

368  PRINT  Q(I) 

369  NEXT  I 

370  READil:  18 

371  FOR  lal  TO  18 

372  READil:  D(I) 

373  PRUT  DU) 

374  NEXT  I 


100’ 

110' 

120' 

130* 

140* 

150* 

160* 

170' 

190* 

-  200' 
220* 
230' 

w  240* 

250* 
260* 
270* 
280* 
290* 
320* 
330* 
340* 

I.  350* 

360' 


I 

i 
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INTEGX  (con tinned) 


575  I8-I8-1 

380  DIN  X(50)  ,G(50)  ,L(  3  *50)  ,S(  1000 ,6) 

390  PRUT 

400  FOR  I«1  TO  72 

410  PRUT  "a"; 

420  NEXT  I 
430  PRUT 

435  PRINT" PLOT  FILE  NAME" 

436  UPUT  T$ 

437  F1LEM2:  T$ 

438  SCRATCH#2 
440  PRUT 
450  GO  TO  990 

460  IF  IfcO  THEN  530 
470  PRINT  "PLOT  SYMBOL"; 

480  FOR  1*1  TO  M 

490  PRINT  TAB(  15*1+3)  ;E$<  I); 

500  NEXT  I 
510  PRINT 
520  PRINT 
530  PRINT  "T"; 

540  FOR  I«1  TO  4 

550  IF  N(I)aO  THEN  600 

560  PRINT  TAB(  15*1)  ;"X("  ;N( I)  ; 

570  Z(I)sN(I) 

580  N2-N2+1 
590  NEXT  I 
600  PRINT 

610  IF  T<T 2  THEN  970 
620  GO  TO  890 
630  GO  SOB  1290 
640  FDR  I»1  TO  N 
650  L(1,I)«H»G(I) 

660  X(I)aX(I)+L(1,I)/2 
670  NEXT  I 
680  T-T+H/2 
690  GO  SOB  1290 
700  FOR  I»1  TO  N 
710  L(2,I)*H*G(I) 

720  X(I).X(I)-L(1fI)/2+L(2,I)/2 

730  NEXT  I 

740  GO  SUB  1290 

750  FOR  1*1  TO  N 

760  L(3,I)«H*GU) 

770  X(I)»X<I)-L<2tI)/2+L(3'I> 

780  NEXT  I 
790  T»T+H/2 
800  GO  SUB  1290 
810  FOB  lei  TO  N 

820  X(I)aX(I)*L(3»I)+<L(1fX)+2*L(2fX)+2*L(3'i)+H*G(X))/6 
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INTEGX  (oontlnuad) 


830  NEXT  1 

WO  T«IMT(  T» 1000004-  .5)/ 100000 
850  H2aH2+1 

860  IF  H2<H1-.5  THEN  970 
870  IF  T<T2  THEN  970 
880  H2>0 
890  JaJ+1 
900  S(J,1)bT 
910»  PRINT  T; 

920  FOR  lal  TO  N2 

930*  PRINT  TAB<15*I>;X<N<I»; 

940  S(J,U1)>X(N(I)) 

950  NEXT  I 

955  PRINT#2:T;X(IW1) 

960*  PRINT 

970  IF  T>T1  THEN  1090 
980  00  TO  630 

990  READ  H,T2,T1,M 

991  NbI8 
1010  READ  HI 
1020  FOR  lal  TO  4 
1030  READ  N(I) 

1040  NEXT  1 

1045  N(1)sI8»1 
1050  FOR  lal  TO  M 
1060  READ  I(I)»E$(I) 

1070  NEXT  I 
1080  GO  TO  460 
1090  PRINT 

1095  PRMT#2:"  1E37  ,  1E37  « 

1100  IF  HbO  THEN  1230 
1110  FOR  lal  TO  15*(N2+1) 

1120  PRINT  •*-■; 

1130  NEXT  I 
1140  PRINT 
1150  PRINT 

1160  PRINT  "VARIABLE*' ,  "SYMBOL" 

1170  PRINT 
1180  FOR  lal  TO  M 

1190  PRINT  **X(";2(I);")">TAB(17);E»(I> 

1200  NEXT  I 
1210  PRINT 
1230  PRINT 
1240  FOR  lal  TO  72 
1250  PRINT  "a**; 

1260  NEXT  I 
1270  PRINT 
1280  GO  TO  1370 

1290  REM  •»•»«  SUBROUTINE  TO  EVALUATE  THE  GUVS 
1292  FOR  lal  TO  18-1 


n  fotM  tsto  r»iN:fo  in  u.s.a. 


4- 


w  UTEGX  (continued) 

1293  GU>  >  X(I+1) 
w  1294  NEXT  I 

1296  G(18)  •  1. 

1298  FOR  1«1  70  18 

v.  1299  G(I8)aG(Z8)  -  D(I)»X(I) 

1302  NEXT  I 
1305  X(I8+1)*0* 

«  1306  FOR  I*  1  TO  17 

1307  X(I8*1).X(I8*1)  «•  Q(I)*X(I) 

1309  NEXT  I 

1320  RETURN 

1330  DATA  .001,0,20,0 

1352  DATA  50,1,0,0,0 

1360  DATA  1.''»,,,2,"#«,3,%,,,4,"«“ 

1370  E NO 


w 
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APPENDIX  D.  COMPUTER  PROGRAM  FOR  THE  TIME 

SOLUTION  OF  THE  SYSTEM  EQUATIONS. 

The  program  on  the  next  page  solves  the  set  of  equations 

x  ■  Ax  ♦  B6 

en  -  C,x 
o  4— 

z  -  Fz  ♦  B0  +  Ke„ 

—  —  —  o 

A 

e  *  Mz 

The  equations  included  are  those  of  the  ideal  observer  and 


rate  sensor. 


1 


IMTEG 


I  , 


V  J 


\ 


100 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
23 0 
240 
250 
260 


1NTEG 

DESCRIPTION — Solves  a  set  of  N  first-order  differential  equations 
of  the  form: 


dX(l)  «  G(I)  f or  I  -  1,2, 


when  the  initial  oonditons  are  known,  by  the  f birth-order  Runge  Kutta 
technique  and  provides  a  plot  of  T  vs  X(N1),  X(N2),  X(N3),  and 
X(N4)  in  that  order  depending  on  selection  in  DATA  statement. 

INSTRUCTIONS— The  functions,  G(I),  are  to  be  inserted  after  line  1350 
and  before  line  300  in  the  following  form: 

1350  G(1)  *  WHATEVER  IT  IS 

1360  G(2)  >  WHATEVER  IT  IS 


270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370' 

380* 

390' 

400' 

410' 

420' 

430' 

440' 

450' 

460' 

470' 

480 

490 

500 

510 

520 

530 

531 

532 

533 

534 

535 
550 


The  entire  mathematical  model  should  go  in  this  section  including 
algebraic  constraints. 

1520  DATA  N,H»T2,T1 ,M 

1530  DATA  X(1),X(2), . ,X(N) 

1540  DATA  H1,N1tN2,N3,N4 

Where  N  *  Number  of  dependent  variables 
H  *  Increment  size 
T2  a  Initial  value  of  T 
T1  *  Final  value  of  T 
M  a  Number  of  variables  to  be  plotted 

HI  a  No.  of  increments  of  H  between  print  outs  (must  be  an 
integer.) 

N1,N2,N3»N4  identify  the  dependent  variables  to  be  printed  out. 

In  line  1530  X(1)...J((N)  a  Initial  values  of  X(1)...X(N).  In 
general,  X(1) ...*X(N)  a  Dependent  variables. 

T  a  Independent  variable 

G(I)  a  Known  functions  a  dX(I)/dT 

I 

•  »  •  »  MAIN  PROGRAM  •  »  «  * 

LIBRARY  "L  .ES***:PL0T1" 

DIM  X(50),G(50),L<3,50),S<500-,6> 

FOR  1*1  TO  72 
PRINT 
NEXT  I 

PRINT" PLOT  FILE  NAME" 

INPUT  T$ 

FILE#2:  T$ 

SCRATCH# 2 
PRINT 
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-2- 


INTEG  (continued) 


4 


w 


f 


I 

t 


V* 


i 

] 


l 


560  GO  TO  1100 
570  IF  M.0  THEN  640 
580  PRINT  “PLOT  SYMBOL"; 

590  FOR  1*1  TO  M 

600  PRINT  TAB<  15#I+3>  ;E$(I) ; 

610  NEXT  I 
620  PRINT 
630  PRINT 
640  PRINT  "T“; 

650  FOR  1.1  TO  4 

660  IF  N(I)»0  THEN  710 

670  PRINT  TAB(15*I)»*X(";N(I);")“; 

680  Z<I).N(I) 

690  H2.N2+1 
700  NEXT  I 
710  PRINT 

720  IF  T<T2  THEN  1080 
730  GO  TO  1000 
740  GO  SUB  1360 
750  FOR  1*1  TO  N 
760  L(1,I)»H*G(I) 

770  X( I)*X( I)+L( 1 ,  I)/2 
780  NEXT  I 
790  T-T+H/2 
800  GO  SUB  1360 
810  FOR  1.1  TO  N 
820  L(2,I).H*G(I) 

830  X(I).X(I)-L( 1 ,I)/2+L(2tI)/2 

840  NEXT  I 

850  GO  SUB  1360 

860  FOR  I«1  TO  N 

870  L(3,I)*H*G(I) 

880  X(I).X(I)-L(2,I)/2*L(3,I) 

890  NEXT  I 
900  T.T+H/2 
910  GO  SUB  1360 
920  FOR  1.1  TO  N 

930  X(I)«X(I)-L(3»I)+(L(1»I)+2*L(2,I)+2BL(3tI)+H*G(I))/6 
940  NEXT  I 

950  T» INK T* 100000+  .5)/ 100000 
960  H2.H2+1 

970  IF  H2<H1- .5  THEN  1080 
980  IF  T<T2  THEN  1080 
990  H2.0 
1000  J.J+1 
1010  S(J,1)»T 
1020  PRINT  T; 

1030  FOR  I>1  TO  N 2 

1040  PRINT  TAB(  15*1)  ;X(N(D) ; 

1050  S(J,I+1>»X<Nd)) 
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1NTEG  (continued) 


1060  NEXT  I 

1065  PRINT#2:T;X( 10) 

1070  PRINT 

1080  IF  T>T1  THEN  1200 
1090  00  TO  740 
1100  READ  N,H,T2,T1 ,M 
1110  HAT  READ  X(N) 

1120  READ  HI 
1130  FOR  1*1  11)  4 
1140  READ  N(I) 

1150  NEXT  I 

1160  FOR  1*1  TO  H 

1170  READ  I(I),£$(I) 

1180  NEXT  I 
1190  GO  TO  570 
1200  IF  H*0  THEN  1320 
1210  FOR  1*1  TO  15*<N2*1) 

1220  PRINT 
1230  NEXT  I 
1240  PRINT 

1250  PRINT  "VARIABLE", "SYMBOL" 

1260  PRINT 
1270  FOR  1*1  TO  M 

1280  PRINT  "X(" ;Z(I) ;")",TAB(  17) ;E$( I) 

1290  NEXT  1 
1300  PRINT 
1320  FOR  1*1  TO  72 
1330  PRINT 
1340  NEXT  I 
1350  GO  TO  1560 
1360  K1* .0253 
1370  K1* .0253 
1380  K2» .04112 
1390  K3* .984566 
1400  K4.-167.3 
1410  G(1)*X(2) 

1420  G(2)*X(3) 

1430  G(3)*X(4) 

1440  G( 4) *-1023 «4*X( 0-4738 .4»X<2)-15653*X<3M75,3*X(4)+1170»X(9) 

1450  G(5)*K1*(X(4)-X(8))fX(6) 

1460  G(6)*K2»<X<4)-X<8))*X<7) 

1470  G(7)*K3t(X(4)-X(8))^X(8) 

1480  G(3)*K4»(X(4)-X(8))-1023.4«X(5)-4738.4»X(6)-156534X(7)-175.3*X(8)+1170«X(9) 
1490  G(9)*1 

1500  X( 10)*13 .33*X(8)+4 .04«X(7)+ .87*X(6) 

1510  RETURN 

1520  DATA  10, .005,0,20,0 
1530  DATA  0,0, 0,0, 0,0, 0,0, 0,0 
1540  DATA  10,10,8,0,0 
1550  DATA  1,"«\2,"#",3,,V,,4,"r 
END 
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